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SUMMARY 

The paper describes a method for guiding a dynamic system through a given set of points. The 
paradigm is a fully automatic aircraft subject to air traffic control (ATC). The ATC provides a sequence 
of waypoints through which the aircraft trajectory must pass. The waypoints typically specify time, 
position, and velocity. The guidance problem is to synthesize a system state trajectory that satisfies both 
the ATC and aircraft constraints. Complications arise because the controlled process is multidimensional, 
multiaxis, nonlinear, highly coupled, and the state space is not flat. In addition, there is a multitude 
of operating modes, which may number in the hundreds. Each such mode defines a distinct state 
space model of the process by specifying the state space coordinatization, the partition of the controls 
into active controls and configuration controls, and the output map. Furthermore, mode transitions are 
required to be smooth. 

The proposed guidance algorithm is based on the inversion of the pure feedback approximation, 
followed by correction for the effects of zero dynamics. The paper describes the structure and major 
modules of the algorithm, and the performance is illustrated by several example aircraft maneuvers. 



1 INTRODUCTION 


Development of a methodology for the design and verification of aviation operations systems is a 
NASA goal. (See ref. 1.) For the purposes of the present report, several interacting components of such 
a system may be identified as shown in figure 1 . They include air traffic control (ATC) and a varying 
number of flights, each consisting of a pilot crew, the aircraft, and the flight vehicle management system 
(VMS). Finally, the ever-present environmental conditions, including weather and hardware failures, 
affect all the other components. 



Figure 1. Major components of an aviation operations system. 


Consider the vehicle management system; its purpose is to provide interfaces between the human 
pilot, the ATC, and the aircraft. The major subsystems of VMS and their interactions are shown in 
figure 2. (See ref. 2.) The human pilot interacts with the VMS at three levels. At the highest level, 
which is knowledge-based (termed strategic in the upper left-hand comer of the diagram), the pilot in 
command generates the flight plan and proposes it to ATC, accounting for weather, route, fuel reserves, 
and other parameters. ATC then modifies this proposed flight plan as necessary to provide separation 
from other traffic, and issues a clearance to the flight. On board the aircraft, the clearance is checked 
for legality and feasibility, and either accepted or renegotiated with ATC, as necessary. The clearance 
is then translated into an executable reference flightpath, and appropriate modes are selected to enable 
the aircraft to fly the accepted clearance. 
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In addition to operational modes such as takeoff, climb, cruise, approach, and landing that are 
selected by the pilot at the strategic level, lower-level modes are selected by the pilot at the tactical 
level. These tactical modes include guidance, configuration, flight control, navigation (state estimation), 
and display modes. 

At the skill level, the pilot is manually controlling the aircraft with little assistance from the 
automatics. 

It may be noted that there is a mixture of both computation and message types of variables. Some 
computations, particularly at the high levels, are discrete or sentential in that they are based on a formal 
logic and they manipulate sentences in that logic, while other computations are arithmetic or numeric 
in the sense that they employ mathematical methods that manipulate real variables. Similarly, the 
communication between subsystems may be by means of both sentences and numbers. In the diagram, 
dashed lines represent flows of sentential or boolean variables; solid lines represent flows of numbers. 
For the present discussion, continuous variables and sampled data need not be distinguished between, 
but both can be referred to as analog variables. Systems with mixed variables, both discrete and analog, 
are referred to as hybrid. It is apparent that a VMS is a hybrid system. 

In the diagram, the aircraft represents the airframe, including engines, servo actuators for all con- 
trols, and all sensors. In this view, the aircraft is the main source of analog variables. It is modeled 
mathematically with differential equations, which are multidimensional, multiaxis, highly coupled, non- 
linear, and time varying. 
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The interaction between the sentential and analog components of a VMS is further illustrated by the 
diagram in figure 3. In this model, the VMS communicates with ATC by means of sentences and with 
the controlled aircraft by means of analog signals. The virtue of the diagram is that it shows that there 
must be an interface, called herein a dynamics translator, which translates sentential commands into 
functions of time, and conversely, analog measurements into sentential descriptions of system behavior. 
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Figure 3. Interface between sentential and analog layers. 


The multiaxis servo is the analog layer. The arithmetic controller closes the loops around aircraft 
disturbances and modeling errors, thereby making the analog layer robust, and it provides accurate 
tracking of input signals. 

Another representation of the aviation operations system is shown in figure 4. For simplicity, only 
the fully automatic mode of operation, in which the flight is controlled completely by ATC and flight 
computers, is shown. A hierarchy of subsystems or layers is in effect, namely, (top to bottom) ATC, 
flight, VMS, path generator, guidance, and regulator. The hierarchy operates in terms of the horizon 
width, persistence, detail, and type. In planning a sequence of way points, ATC has the widest horizon, 
the plan changes infrequently (minutes to hours), it is very coarse, and the planning algorithm is of the 
discrete type. The path generator, in converting the waypoints into control points (described later in the 
report) has a narrower horizon than does ATC; replanning occurs more frequently (seconds to minutes); 
there is a moderate amount of detail in the plan; and planning employs a mix of discrete and continuous 
methods. The guidance, which transforms the path into the corresponding reference motion of aircraft 
state and control ( x 9 ,u 9 ), has a narrow horizon with only several seconds of “look-ahead” capability; 
replanning occurs frequently (several seconds), and the planning employs mainly continuous methods. 
Finally, at the bottom of the hierarchy, the regulator transforms the tracking error x e = x p — x 9 , which 
is the difference between the estimated plant state and the guidance state, into corrective control that is 
combined with the guidance control u 9 to produce the plant control vP. (The negative of the standard 
servo error is chosen so that x 9 may be treated as the origin.) The regulator has the narrowest horizon, 
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Figure 4. Multilayer model of aviation operations system. 


without any look-ahead capability; the shortest persistence (milliseconds); and the greatest amount of 
detail. The algorithm is based on continuous (and sampled-data) methods. 

The environment injects disturbances into the aviation operations system at many levels. The 
standard procedure for reducing the effects of these disturbances is to employ feedback. Centralized 
feedback in which the hierarchy is ignored is always possible, but it leads, in all but the simplest cases, 
to very complex and largely unsolvable performance analysis problems. It seems that a more promising 
approach is to structure the feedback hierarchically and to match at each level the variable type of the 
feedback with the algorithm type of the level. If the feedback can be layered in this way, then the 
complete system will consist of a hierarchy of closed-loop subsystems. Each subsystem has the servo 
structure in the sense that its function is to closely follow its command in spite of disturbances. In 
such a completely layered structure, the type of command and disturbance depends on the type of the 
layer. The expected analytical advantage of such a structure is that each servo has a relatively simple 
input/output description, and communication between servos is simpler. This is not passive modeling — 
it is active design to make it so. Part of the freedom available in the design process is used to structure 
the system; the remainder is used for control. The present working hypothesis is that such complete 
layering is possible for the aviation operations system. Consider the block diagram shown in figure 5. 
In the figure, the ATC waypoints are represented by a waypoint card sequence (WCS) that describes the 
waypoints. The discrete event control system operates on the waypoint card sequence, the pilot input, 
and the system behavior represented by the status card sequence (SCS), and produces the control card 
sequence (CCS). The control cards are the only means for controlling the rest of the system below this 
level. That is, at this level the CCS is the control variable. Furthermore, the SCS is the only means for 
observing the behavior of the system below this level. 
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Figure 5. Closed-loop hybrid flight control model. 


The guidance output trajectory generator translates the control cards into appropriate continuous 
functions of time that define the commanded motion of system output y c , which may include, for 
example, the runway coordinates of aircraft position vector, the position of flaps, and other variables. 
This generator also specifies the operation mode //, which specifies some of the details (such as control 
mode) of just how the maneuver given by y c is to be executed. 

The guidance state trajectory generator expands y c into a complete state trajectory and correspond- 
ing control (x g ,u g ) that satisfies the state equation of the aircraft model and the output map. The 
aircraft is then regulated to that reference state. Finally, the performance monitor reports any deviations 
from the expected motion by means of the SCS. 

This report is concerned only with the transformation of CCS into the corresponding guidance 
( x g ,u g ). The guidance algorithms that are developed are quite general; they apply to both civilian 
and military operations, including remotely piloted vehicles, and to a variety of aircraft types such as 
conventional takeoff and landing (CTOL), short takeoff and landing (STOL), vertical takeoff and landing 
(VTOL), helicopters, and tilt rotors. 

The guidance problem considered in the present report may be described as follows. The control 
point table describes the control points. For example, figure 6 shows a 90-degree turn that is composed 
of several segments, which are specified by control points. The straight west-to-east line segment 
terminates on the control point c The next segment is a variable curvature transition to a circle at 
c k+l- The circular arc, terminating at cj. + 2 , is followed by a variable curvature transition terminating 
on a north/south line at c/ (+ 3 . The CCS provides the x,y coordinates at each control point. In addition, 
it may specify the duration or the arc length of the segment, as well as a set of consistent terminal 
conditions for some of the many additional variables such as altitude, climb rate, ground speed, airspeed, 
heading, thrust level, and flap deployment. The CCS also provides higher level instructions such as 
the type of segment (linear, circular, or helical), and it specifies the operating mode for the segment. 
The operating mode provides such high-level instructions as (1) which of several possible state space 
representations is to be used; ( 2 ) which of several possible state-to-output maps is to be used; and 
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Figure 6. Typical maneuver specified by means of control points. 

(3) which of several control modes is to be used. The operating mode is discussed in detail in the 
main body of the present report. At this point we merely note that there may easily be over a hundred 
possible operating modes. 

Each such mode (denoted by p) defines a different state space model of the same underlying 
physical process. That is, each operating mode defines a particular state space (X M ), control space (U M ), 
parameter space (P,J, output space (Y /L ), state flow (/^), and output map (h^): 

t e R, x e X /x , u eU^, p e Pfi, y e Y fl 

X = fn(x,u,p) ( 1 ) 

y = M x ) 

The dimension of the state space X^ remains constant; the meaning of the coordinates is mode- 
dependent. All controls are distinguished into two types: The controls used by the regulator are 
called the active controls (or controls, for short) u\ the rest of the controls are called parameters p. The 
dimensions of the active controls and output spaces, U ^ and Y^, respectively, while always equal to 
each other, may be mode-dependent. For example, in one control mode the throttle may be an active 
control that is used by the control system to execute a particular maneuver; in that case the throttle is 
included in U^. On the other hand, occasions arise where the throttle is not available for active control, 
as when set to idle. In that case the throttle is considered to be a parameter coordinate in P ^ and not a 
control coordinate in U fi . Similarly, both airspeed and climb rate may be controlled outputs, in which 
case both coordinates are in Y^. But there is a mode in which the climb rate is not controlled explicitly; 
then it is not one of the coordinates in Y^. Some variables, such as the wind, are always in P^. 

In general, each such model M fl is multidimensional, multiaxis, cross-coupled, nonlinear, and, 
because of p, time-varying. The problem is to find for every segment k the evolution of state and 
control, 

(x 9 ,u 9 )(t), t € I k (2) 

such that they satisfy both the state and output equations. That is, for the duration of the segment 
t € ifc, and for the given desired output y 9 , 

x g {t) = f^[x 9 (t),u 9 (t),p{t)} 3 

y 9 {t ) = h^[x 9 (t) ] 
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Furthermore, the segments must be patched smoothly together. A solution of this class of problems is 
the subject of the present report. 

The report consists of three major sections. First, a general model of an aircraft is developed, 
and the typical structure of the force and moment generators is described. The Harrier VTOL fighter 
aircraft is used as an example. Next, the command generator in which the aircraft maneuvers are 
defined is developed. The maneuvers are represented by a sequence of control points that are linked 
with polynomials. Several useful representations of aircraft attitude are also presented. Finally, the 
multimode inverse algorithm, which generates the initial state and the control time history that produce 
the desired maneuver, is developed. The methodology is illustrated by means of several numerical 
examples. 

The aircraft model used in the report is quite realistic. As a result, it is relatively complex. The 
state space is high-dimensional and, because of rotations, not Euclidean. The model is multiaxis (up 
to six), cross-coupled, and nonlinear. Furthermore, there is a need to employ several representations, 
output maps, and control strategies. The methodology of dynamic forms developed in “Dynamic Forms 
Part I: Functions” (ref. 3) provides effective means for the management of this complexity. A dynamic 
form of order k of a variable of time is simply the variable together with its first k time derivatives. 
The essential aspects of the methodology are reviewed and employed throughout the present report. 
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2 AIRCRAFT MODEL 


In this section a state equation describing the open-loop dynamic behavior of an aircraft is developed 
in several steps. This equation forms the basis for a variety of state space models that are developed 
and used in the remainder of the report. 


2.1 Rigid Body 


The core of the model is a rigid body moving in three-dimensional, inertial space, subject to 
the gravitational, aerodynamic, and propulsion forces. In this report the dynamic effects due to Earth 
rotation and nonflatness are ignored. This approximation causes less than 0.004 g acceleration error 
for speeds up to 600 mph, and such errors are consistent with the assumed modeling accuracy of the 
aerodynamic and propulsion forces. Similarly, the aircraft mass, location of the center of mass, and the 
moments of inertia are assumed to be quasistatic in the sense that their time derivatives do not appear 
in the model. Furthermore, the effects of rotating parts of the aircraft are ignored. All such translational 
and angular acceleration errors that are of the order of 0.004 g and 0.002 radians/sec 2 , respectively, are 
ignored since their effects can be easily controlled by means of a realistic regulator. 

Let V be a local (runway) right-handed, orthonormal axis system in which 7% points in the^direction 
of the local vertical (down), f\ points in the direction of the runway, and r*2 = f 3 x rj . The r system 
is attached to the runway. The runway coordinates of vectors are denoted by the subscript r. 


Several variables are fundamental to the following discussion. The position, velocity, and accel- 
eration of the aircraft center of mass relative to the runway are denoted by f, v, and a, respectively, 
the average wind velocity at r by w, and the relative air velocity by v 0, = v — w. The aircraft mass 
is denoted by m, the total nongravitational force acting on the aircraft by F, the specific force by 
/ = F/m, and the gravitational force by mgr , 3 with scalar, constant g. Since r is assumed to be 
inertial, the standard form of the translational equations of motion follows: 

f r = v r (4) 

V r = a r = fr + g8 3 

where for i = 1, 2, 3 is a column matrix with zero entries except in row i, where the entry is 1. 

4-4 

Another frequently used axis system is attached to the aircraft. In this ^body) system, b, b\ 
points tail to nose, 63 is in the aircraft plane of symmetry and perpendicular to b\, pointing down, and 

b 2 = 63 x b-[ points along the right wing. The b system is attached to the aircraft center of mass. 


The attitude of the aircraft relative to the runway is defined by the direction cosine matrix C^ r so 
that, for example, the body coordinates of the aircraft velocity v are given by 

Vf) = Cfo r v r (5) 

The inverse of the direction cosine matrix equals its transpose: C^. 1 = Cf r = C^. The body coordinates 

of the angular velocity of b relative to r are denoted by u)^, so that 

Cbr = SMC br (6) 
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where, for any three-dimensional column matrix x. 


(7) 


S(x) = 

The function S represents the vector cross 
x r = —S(y r )z r . 


( 0 2:3 -X2 \ 

-xz 0 x\ 
x 2 ~ x \ 0 ) 

product. For example, x = y x z is represented in V as 


The angular momentum and the total torque about the aircraft center of mass are denoted by L and 
N, respectively. Since r is assumed to be inertial, 

L r = N r (8) 

where the overdot denotes the time derivative of the (scalar) elements of L r . But Lb = Q, r L r , so that 
after applying equation (6), 

Lb = Nb + s iNbrb)Lb (9) 

The body coordinates of the moment-of-inertia operator about the center of mass are denoted by Jf t so 
that, ignoring rotating components, 

Lb — Jb^brb 00 ) 

and, since the time derivatives of Jb are ignored, the Euler equations (refs. 5 and 6) become 


Jb^brb = N b + S{u)brb)Jb^brb 

or, equivalently, 

O brb n b "b '^b ^ i^brb) ^b^brb 

where the specific torque rib — For convenience, let 

m b = ^b "b ^ b (idbrb) Jb^brb 

The parameter mb is the specific moment. 

In summary, the equations of motion of a rigid body in inertial space follow: 

V f = -- Vf 

v r = C r bfb + gh 

Cbr = S (oJbrb)L'br 
O brb = m b 

Next a model of the specific force fb and moment mb is developed. 


( 11 ) 

( 12 ) 


03) 


(14) 


2.2 Force and Moment Generators 

In this section the structure of the function relating the aerodynamic, propulsion, and reaction 
control forces and moments to other process variables is developed. The special case of the Harrier 
aircraft is considered, but the formalism applies to other aircraft types. 
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The YAV-8B Harrier (ref. 4) is a powered-lift vertical/short takeoff and landing (V/STOL) transonic, 
light, attack aircraft. The engine is a single turbofan, whose thrust may be diverted by means of a 
nozzle system. The directed thrust supports the aircraft in low-speed flight. Aerodynamic lift supports 
the aircraft in high-speed, conventional flight. The general features are shown in figure 7. The engine 
thrust vector is controllable in the longitudinal plane by means of the throttle, which controls the 
magnitude, and the nozzle, which controls the direction. The control of aircraft rotation depends on 
the flight condition. At high speed, the torque vector is generated aerodynamically by means of the 
conventional ailerons, elevator, and rudder. At low speed, when the aerodynamics is ineffective, the 
torque vector is generated by means of a reaction control system (RCS). In the RCS, a portion of the 
engine gas is ducted to the extremities of the aircraft, where the directed flow produces reaction forces 
whose magnitude is controlled by means of valves. The valves are connected to the aircraft roll, pitch, 
and yaw controls. In the transition regime, there is a smooth blend of aerodynamics, powered lift, and 
RCS. Note that, since a request for moment may produce unwanted forces as a side effect, there may 
be zero dynamics. 



The algorithm YAVSB adapted from reference 4 for computing the total force and moment is 
shown in figure 8. At the top level, the algorithm implements the function 

YAVSB : R 19 -► B? x R 3 (15) 

The input consists of moment controls u m G .R 3 , power controls u p G i? 3 , configuration controls 
6 /? 3 , air velocity v% G i? 3 , angular velocity € -R 3 , pitch angle 0, altitude h, air temperature 
and aircraft weight W. The output consists of the specific force /& G i? 3 and the specific moment 
mf, G .R 3 . In this case of the Harrier, the power control u p consists of the throttle, empty slot, and 
nozzle. In general, the middle term could be a side-force controller. The configuration control uj 
includes the flaps, aileron droop, and the speed brake. 
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Figure 8. Force and moment generation algorithm for the YAV-8B Harrier aircraft. 


At the next level down are seven major subalgorithms. Secondary variables are computed first 
from the input. The aircraft mass, the location of the aircraft center of mass, and the moment of inertia 
matrix are computed in MASS DAT A. The air density, pressure, and speed of sound are the outputs of 
ATMOS, which models the atmosphere. The airspeed, angles of attack and sideslip, dynamic pressure, 
and Mach number are provided by AIRDATA. The preliminary computations are completed in the 
FANOUT algorithm, in which the controls ( u m ,u p ,Uf ) are interpreted in terms of the engine throttle, 
the nozzle angle, the eight reaction control valve openings, and the seven aerodynamic control surfaces. 

The details of the actual physics of the process are represented by the three large algorithms AERO, 
ENGINE, and RCS, which compute the aerodynamic, engine, and reaction control system torque and 
force. In the case of Y AV8B, there is a one-way interaction between engine and aerodynamics through 
the jet effects and a two-way interaction between the engine and the RCS due to bleed gas ejected 
through the reaction control valves. The pitch angle 9 and altitude h are used in AERO to compute 
the ground effect. Other aircraft types, such as other VTOL types, powered-lift STOLs, helicopters, 
tilt rotors, and CTOL aircraft have similar algorithms describing the corresponding force and moment 
generators. The action of the algorithm is denoted by 

Vb = F ( z b ) (16) 

where the input and output are given by, respectively, 

z b = ( u m, Up, Uf, v%,u brb , 9, h, W,9 T ) e (-R 3 ) 5 X (i?) 4 

Vb = (fb,ui b )e(R 3 ) 2 

Henceforth it is assumed that the algorithm is given for the particular aircraft being considered. 
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In general, the body coordinates v b of the air velocity, which is an entry in the input z b , may be 
computed as follows 

v b =c br (v r -w r ) (18) 

where w r € R 3 represents the runway coordinates of wind. Another entry, the altitude h, is the negative 
of the third coordinate of position: 

h = — r,- 3 (19) 

and the pitch angle 6 is the middle Euler angle in the following expansion of the direction cosine matrix 

C br = E l (4>)E 2 (6)E^) (20) 

Therefore, z b may be computed from (u,x w ,p), say 

z b = Z(u,x u ,p) (21) 

where u, x^, p are the combined control, rigid body state, and parameters, respectively, defined as 
follows: 

u = ( u m , Up , uj ) £ U C (R 3 ) 3 

Xa; = ( r r , Vr, c br , u brb ) € Xu C ( R 3 ) 2 x 50(3) X R 3 (22) 

p = (w r , W, 0 T ) e P C R 3 x R x R 

The combined algorithm is denoted by F b = ( f b ,f b n )■ Consequently, equations (14) may be 
written as the rigid-body state equation, as follows: 


V f — Vf 

v r = C rb fl (u, x-uj , p) + g6 3 
C b r = S{^brb)Cbr 

^brb f b (UfXuiP') 

with the arguments defined by equations (22). 


(23) 


2.3 Simulation Model 

In general, the nine controls ( u m ,up,ut ) are not directly accessible. Rather, intervening dynamics, 
including position and rate limits, exist between the accessible servo inputs (u£ t , u £ , itj) and the outputs 
(u m ,u p ,uf). However, often the control servos may be adequately modeled by the usual second-order 
systems of the form 

— = — cy ^ 7 (24) 

u c s 2 4 - 2(u n s + 

as shown figure 9. The servo controls (u„ v Up, iPj) are connected to the inputs (u£ t ,u£,Uy-) by the 
usually decoupled servo control laws: 

u^ = a m (u m ,u^ n ,uf n ) 

Up cr p (u p , Up, Up) (25) 

Uf = a f(uf,iif,u C f) 
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such that there is good tracking for slowly varying inputs that remain within their position limits: 

(u m ,u p ,Uf) ta (u^u^uj) ( 26 ) 

The other inputs to the model are the parameters, namely, runway coordinates of wind w r , aircraft 
weight W, ambient temperature Qj, and the rigid-body state x u . The block Xf> produces the variable 

x h, which is the input to the force f[ and moment algorithm. The body coordinates of force f\ } are 
transformed into the runway coordinates f r and the effects of gravity are added, resulting in the aircraft 
acceleration a r , which is integrated to produce aircraft velocity v r , and another integration produces 
aircraft position r r . 

In the rotational part, the body coordinates of the angular acceleration of the aircraft about its 
center of mass, = m^, are integrated to yield the body coordinates of angular velocity relative to 
the runway, a and another integration according to equation (6) results in the direction cosine matrix 

Ci r , locating the body axes b relative to the runway axes r . 

It should be noted that the presence of the direction cosine matrix Q> r causes the state space to be 
nonflat, but a small angle patch of the state space is a subset of Rr Q . The analysis of the large angle 
behavior of the system is complicated by this relatively large dimension and type of the state space. 
The analysis is further complicated by the nature of the force and moment generator, which is generally 
high-dimensional, multiaxis, cross-coupled, and nonlinear. 
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In a typical real-time, piloted simulation, the algorithms shown in figure 9 are connected to several 
other systems, such as a cab that provides to the pilot translational and rotational motion cues, and a 
cockpit that provides the pilot access to the controls and visual, tactile, audio, and electronic displays. 
The simulations are quite realistic. In fact, in many cases the control software is loaded into the 
simulation computers, where it is checked out by the test pilot and then loaded without changes and in 
exactly the same way into the flight computer and flight tested. 

The remainder of the report assumes a high-fidelity force and moment algorithm = (f[, f ™). 
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3 OUTPUT TRAJECTORY COMMAND GENERATOR 


The purpose of the command generator is to provide the interface between the discrete and analog 
levels. As already noted, the signals in the discrete level are best described as sentences, which compose 
the current plan of action generated by the automatic reasoning system. The signals in the analog level 
are real functions of time. At the detailed level these signals are, of course, not analog but sampled 
data. However, since the present variables are relatively slowly varying, while the usual sampling rate is 
high, sampling is largely ignored, and the signals are analog. The high-level plan to be executed by the 
aircraft may be considered to be a series of maneuvers. Furthermore, there is a hierarchy of maneuvers 
in the sense that a particular maneuver may be further parsed into a more-detailed concatenation of 
simpler maneuvers. The maneuvers have both qualitative and quantitative components. In this chapter, 
a method for converting the commands of the sentential (discrete) level to the analog level is developed. 
The interface is realized by means of a sequence of records (control cards) that specify the complete 
maneuver. In the next chapter, these maneuvers are converted into detailed guidance commands, which 
specify the nominal evolution of system state and controls. 

3.1 Aircraft Maneuvers 

The type of maneuvers that the automatic system should be able to execute is illustrated by the 
following examples. 

Maneuver M\ = M x (Ax, At). In this maneuver, the aircraft is required to step forward a given 
distance Ax in a given amount of time At. Furthermore, the aircraft is in hover at both the beginning 
and the end of the maneuver. This relatively high-level maneuver may be parsed into a sequence of 
more elementary submaneuvers, as follows. 

M\ = hover • accelerate • coast • decelerate • hover 

The specific case of stepping forward 300 feet is shown in figure 10, where motion is from right to 
left. In panel (a) time is zero and the aircraft is in hover with zero pitch attitude and with nozzle at 
90 degrees. The control mode is VTOL, in which both the throttle and the nozzle are available for 
longitudinal force control. The next submaneuver is to accelerate to 8.2 ft/sec 2 , holding pitch attitude 
at zero degrees. At (b), 4.5 seconds later, the position, velocity, and acceleration are, respectively, 
16 ft, 15 ft/sec, and 8 ft/sec 2 ; the nozzle angle changed to 74 degrees to provide the horizontal force 
needed for the specified acceleration. The next submaneuver (panel c) is to coast at 30 ft/sec for two 
seconds. Finally, the aircraft is commanded to decelerate to zero velocity. It may be noted that, while 
the motion of the aircraft center of mass is symmetric about the midpoint (c), the control is not, because 
the available nozzle range is from 5 to 100 degrees. Therefore, the control mode has to be switched 
to one in which the nozzle becomes a parameter (in the present example driven to 90 degrees), and 
the control action flows through the pitch axis. This control scheme of tipping the thrust is typical 
in the control of helicopters. For this reason this control mode is denoted as HEL. In panel (e), at 
time 18 sec, the aircraft is hovering over the point 300 ft from the starting point, as specified by the 
macromaneuver M^OO, 18). 
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Figure 10. The step forward maneuver, M x (300, 18). 


Similar maneuvers in the lateral My (Ay, At) and vertical M z (Az, At) directions are, of course, 
also of practical interest, as are simultaneous combinations, M r (Ar, At). 

Maneuver M 2 = Mq(A 0, At). In this maneuver, the aircraft center of mass remains fixed while 
the pitch angle is changed by the amount A9. The duration of the maneuver is Ai. The case of 
Mq( 81,12) is shown in figure 11, where the drawings are offset for clarity. Initially, the aircraft is 
in the VTOL mode and the pitch is zero. During the maneuver the control mode changes to vertical 
attitude takeoff and landing (VATOL), that is, a “tail sitter,” in which the force control flows through 
thrust, pitch, and yaw channels. The nozzle becomes a parameter driven to 9 degrees, and the maneuver 
includes a change of coordinates. The near-horizontal attitudes are usually represented by means of the 
Euler angles in the standard sequence of roll-pitch-yaw. That representation becomes singular in the 
vertical attitude. On the other hand, the yaw-pitch-yaw sequence is nonsingular in the neighborhood of 
the vertical attitude but becomes singular for near-horizontal attitudes. During the maneuver the control 
system must smoothly patch the two representations. 



Figure 11. Transition from VTOL to VATOL maneuver, M 2 . 
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Similar maneuvers may be defined for the other two axes, namely, roll and yaw. Whether or not 
a given maneuver is executable depends on the available control degrees of freedom. If an aircraft, 
unlike the Harrier, does have a direct three-axis force controller, then translation and rotation become 
independently controllable. In any case, a macromaneuver may consist of simultaneous translation and 
rotation submaneuvers. For example, a vertical step together with a large pitch change would look as 
shown in figure 11 . 

Maneuver A /3 = FRISBEE. This only partly whimsical example illustrates the multiaxis, 
multisegment nature of the maneuvers. It is assumed that an accurate force and moment model exists 
in the unusual flight conditions. As shown in figures 12 and 13, the aircraft is taken from steady hover, 
through a climbing transition, to intermediate speed flight. In addition, the aircraft executes a flat, 360- 
degree yaw maneuver. As before, the motion is from right to left in the figures. Initially, in panel (a), 
the Harrier is in steady hover, at which time the command is given to accelerate speed, altitude rate, 
and yaw. Some of the variables at 7 sec into the maneuver are shown in panel (b). For example, the 
yaw and yaw rate are 90 degrees and 30 degrees/sec, respectively; the nozzle is in vertical position, the 
required acceleration being provided by tipping the thrust by means of the roll angle of —6 degrees. In 
panel (c), the aircraft is flying backward, the acceleration being provided by the forward position of the 
nozzle at 96 degrees, and the roll angle is back at zero. The speed and climb rate are, respectively, 26 
and 15 ft/sec. In panel (d), the roll angle is not quite symmetric in comparison with panel (b) because 
of growing aerodynamic effects. By 18 sec, panel (e), the yaw maneuver has been completed; the speed 
is 50 ft/sec, and the nozzle is at 82 degrees. The Harrier is in longitudinal accelerated, climbing flight. 
The transition from VTOL to the CTOL mode is shown in figure 13. Up to this point the control mode 
was VTOL, the pitch angle was held at zero, and the flap was set to a maximum of 61 degrees. At 
30 sec into the maneuver, panel (g), the control is changed from VTOL to CTOL, in which the lift is 
generated mainly by means of the angle of attack. The nozzle is commanded to its minimum, and the 
flaps are commanded to retract to 25 degrees. The flaps have to be partially retracted to satisfy the flap 
placard, which limits the allowable aerodynamic forces on the flaps for structural reasons. Finally, in 
panel (i), the aircraft has transitioned to conventional flight at 215 ft/sec speed and 15 ft/sec climb rate. 
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Figure 12. Low-speed part of the Frisbee maneuver, M 3 . 



Figure 13. Conventional-flight part of the Frisbee maneuver, M 3 . 


Maneuver M 4 = ROLL. This maneuver, shown in figure 14, illustrates the case in which the 
definition of the system output changes during the maneuver. In panel (a) the fighter is shown in 
horizontal, high-speed flight and in the CTOL control mode. The output is the triple of Cartesian 
coordinates of velocity. The aircraft is commanded to acquire a given vertical velocity, and then to 
switch to the output mode VYZ in which the variables being controlled are the speed (v) and the lateral 
(Y) and normal (Z) components of force. Then the command is given to unload the force by driving 
its component perpendicular to the flight path to zero, as in panel (c), at which point the aircraft is in 
near free fall, with thrust being used for speed control. Meanwhile, the maneuver analogous to M\ is 
executed in the roll axis: 


M^ = hold • accelerate • coast • decelerate • hold 

to accomplish the 360-degree roll. Finally, the force is reloaded and level flight is recovered, as shown 
in panel (g). 
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Several characteristics may be extracted from the preceding examples: 

• There appears to be a hierarchy of maneuvers so that a maneuver at one level may be parsed 
into a sequence of simpler maneuvers from the next lower level. 

• In general, the parsing may not be unique. 

• The maneuvers are described both qualitatively and quantitatively. 

• It is useful to characterize the qualitative behavior of the aircraft by means of operation modes, 
such as control modes, coordinates modes, and output modes. 


The parsing problem, which is, as yet, unsolved, is the subject of current research. The present report 
is concerned only with the representation and generation of the elementary maneuvers. 
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3.2 Polynomial Segments 


Consider the example maneuver M\ shown in figure 10; the objective is construction of a time 
history of position y(t) such that the aircraft accomplishes the step forward in the allotted time and 
such that it is flyable, that is, the aircraft motion is consistent with the existing dynamic and saturation 
constraints. Position is, of course, linked to acceleration by a pair of integrators. Usually, acceleration 
is not a directly accessible control variable. For example, in panel (d) of figure 10, acceleration is 
accomplished by means of rotation, which inserts two more integrators in the control path. Angular 
acceleration is essentially directly accessible because the servos controlling the torque generators have 
a high bandwidth. However, since the servos are rate limited, one more integrator is inserted in the 
control path to obtain continuity of the control variable. Similar considerations apply to other controls. 
As a result, in order to be a possible motion of aircraft position, y and four of its derivatives must be 
continuous at all times. Consequently, the special case of a string of integrators shown in figure 15 
with y = zi is a model of the single-axis-motion command generator. The /c™ time derivative of .zq is 



u x 5 x 4 x 3 X 2 xi 


Figure 15. Single string of integrators as model for single-axis position generator. 
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Since £j(0) = /, £2(0) = 0, and E 2 (t) is nonsingular for t > 0, 


c 2 = E 2 1 {T)[x(T)-E l (T)x( 0)] 

Therefore, the polynomial with n continuous derivatives that satisfies the boundary condition [x(0), x(T )] 
is given by 

P[t t T,x(0),x(T)] = [E\(t) - E 2 (t)E^(T)E 1 (T) ]x(0) + E 2 (t)E^(T)x(T) (32) 

which is the optimal solution for a string of n integrators with the quadratic penalty on the input: 

J = J jx(” +1 )] 2 (t)cIt (33) 

The polynomial solution has the following two computationally desirable properties: Suppose that it is 
necessary to insert a new control point 0 < T 1 < T inside the current interval. That will not change 
the solution, since for t>Ti 


P[t,T,x(0),x(T)] = P{t-T h T 2 -T!,P[Ti, T,x(0),x(T)],x(T)} 
Furthermore, as can be easily verified, the time variable scales, so that, 

E\(t/T) = DjE\{t)D^ , E 2 (t/T) = T~ n D T E 2 (t)D^ 1 


where 


Dj = diag 1 , T , . . . , T h 



and, therefore, 

P[t,T,x(Q),x(T)] = Dj l P\t/T, 1 ,D t x(0),D t x(T)] 
Of particular practical interest is the case of a step of size x\ (T): 


x(0) — 0 

x(T) = *i(T)«! 


(34) 

(35) 

(36) 

(37) 

(38) 


where is a column of zeros except in row i, where it is one. The corresponding solution is 

ipi[t,T,xi(T)] = x\(T)D^ 1 P{t/T, l,0,6i) (39) 

Consequently, properties of the function t — * ijj[t, T, x\(T)\ can be obtained from the normalized 
function r — + ^(r, 1, 1). For example, the maximum absolute value of each coordinate for t € [0, T], 
which is of considerable practical interest, is related to that over r € [0, 1] by the simple relation 


max\ipi[t,T, xi(T)}\ = x\(T)T l t max\'il>i(T, 1, 1)| 


(40) 


So, a unit step in position for T = 4, which is shown in figure 16, may serve as a standard. The 
upper panel shows x\ = (solid) and x 2 = ±\ = (dotted); similarly, the middle panel shows 


(x3,X4) 


J 2) _(3) 

x j , 


, and the bottom panel shows (X5 ,u) = 


x 


(4) _(5) 
1 ’ x \ 


The absolute maxima 
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Figure 16. Unit position maneuver for T = 4. 


[l^^lmax* j — 0, . . . , 5j for a unit step in position for any duration T may be computed from the 
corresponding maxima for the case of T = 4 by means of equation (40). The first row Mqj(T ) in 
table 1 gives these maxima for a unit step in position aq that is accomplished in time T. The second row 
in table 1 represents the maxima for the unit velocity step. For this maneuver the boundary conditions 

are given by x^(T) = 1, Xj(T’) = 0, i = 2, ...5 and x\ (T) is free — it is not matched. Similarly, for 
the other rows in the table the state coordinates to the left of 1.000 are not matched; they are “don’t 
cares.” 


Table 1. Maxima M^AT) of elementary maneuvers 



X 

V 

a 

a 

a 


X 

1.000 

2.461T -1 

9.371T -2 

78.74T -3 

622. 5T~ 4 

15120T -5 

V 

0.500T 

1.000 

2.188T -1 

7.511T’ -2 

52.50T -3 

840. OT -4 

a 

0.143T 2 

0.500T 1 

1.000 

1.875T -1 

5.773T -2 

6O.OOT- 3 

a 

0.033T 3 

0.150T 2 

0.500T 1 

1.000 

1.500T -1 

6.000T -2 

a 

0.008T 4 

0.042T 3 

0.167T 2 

0.500T 1 

1.000 

l.ooor ^ 1 


The unit velocity v maneuver for T = 3 is shown in figure 17. The absolute maxima in this case 
are given by the second row (first derivative) M\A3). For example, |a;( 5 )| max = 840/81. The unit 
acceleration a maneuver for T = 2 is shown in figure 18. 
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Figure 17. Unit velocity maneuver for T = 3. 



time (sec) 

Figure 18. Unit acceleration maneuver for T = 2. 


So far the states were assumed to be unconstrained, but in practice, there may be constraints. 
Thus, it may be necessary to stay within a speed limit, v < v max , and impose limits on the acceleration 
and higher derivatives. For example, in the case of the Harrier in the helicopter (HEL) control mode, 
forward acceleration is accomplished by pitching the aircraft forward, so that, approximately, the forward 
acceleration a r \ = - g6 , where g is the acceleration of gravity and 6 is the pitch angle in radians. 
Situations arise when it is desirable to limit the pitch angle and three of its time derivatives to some 
given values. Pitch acceleration is limited by the capabilities of the pitch moment generator, and pitch 
acceleration rate must be consistent with the rate limits, which may be imposed on the pitch channel 
actuator servo to prevent rapid hard-over failures. That means that the acceleration and higher derivatives 
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should be constrained as follows: 



x^\ \ 


( @max \ 


*< 3 >l 

<9 

@max 

i 

x( 4 )| 

@max 

i 

VI 

x( r> )| J 


i d(3) J 

\ Umax / 


For example, typical constraints, including a speed limit for the Harrier in the low-speed HEL mode, 
may be 


/|x (1) |\ 


/30\ 

|x( 2 )| 


12 

|x(^| 

< 

12 

|x(^| 


12 

V |x< 5 )| / 


\ 24 / 


(42) 


These limits would be violated if, for example, the maneuver M\ = step forward in 18 seconds were 
to be executed as a single segment as shown in figure 19. The maneuver is very smooth, with small 

x^ and but the speed v = x^ exceeds its limit of 30 ft/sec. Next an approach for satisfying 
constraints is considered. 



Figure 19. Unconstrained realization of the maneuver M\. 


3.3 Control Cards 

As is well known, for reasonable situations with constraints, the optimal trajectory is a concatenation 
of segments that are either interior or on the constraining boundary. Strict optimality is not pursued 
further since the structure of the overall control system is of primary concern. However, the insight that 
the optimal solution consists of a sequence of control points that are connected by polynomials is used. 
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In the present case of a string of integrators, in order to stay on the (constant) boundary X{ = xi max , 
the higher derivatives of Xj must be zero. Consequently, the unconstrained elementary maneuvers listed 
in table 2 are considered. 


Table 2. Elementary maneuvers 


a^v.T) 

o$(a,T) 

4(d,T) 


get steady position 
get steady velocity 
get steady acceleration 
get steady accel. rate 
get steady accel. accel. 


x in time T 
v in time T 
a in time T 
a in time T 
a in time T 


The lower index on <7 indicates the first derivative that is not to be ignored; the upper index 
indicates the highest derivative that must be matched. Thus, o\ indicates that position and velocity are 
to be ignored (not matched), and all derivatives of acceleration up to order 3 must be zero at the end of 
the interval. Many maneuvers of practical interest may be parsed into the these elementary maneuvers. 
For example, the overall goal of the step forward maneuver M\ may be achieved in several different 
ways: 

(1) The maneuver may be achieved as a single segment (shown in figure 19) that is given by 

Mi = 4(Ax,r) (43) 

( 7 ) 

The maxima are given by the first row in table 1: | max = AxMqj(T), Ax = 300 ,T = 18. 


(2) The maneuver may be achieved as a symmetric concatenation of three velocity segments (shown 
in figure 20 for T\ =8, T 2 = 2) given by 


Mi = a\{v,Ti)<j\(v,T 2 )o-\{Q,Ti) 

(44) 

where the duration of the maneuver is 


T = 2T\ 4- T 2 

(45) 

and size of the position step is the product 


x = (Ti +T 2 )v 

(46) 


The maxima of velocity and higher derivatives for this maneuver are given by the second row in table 1, 

\vW\max = vM lM i(Ti) (47) 


(3) The maneuver may also be synthesized as a symmetric concatenation of seven elementary 
acceleration segments (shown in figure 21 for T\ =4, T2 = 0, T3 = 2) as follows: 

M 1 = 4(a,T 1 )4(a,r 2 )4(0 ) T 1 ). ( 7f(0,T3).ai(-a,ri) ( 7f(-a ) T 2 )^(0,T 1 ) (48) 
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The duration of the complete maneuver is 



T = 47\ + 2 T 2 + r 3 

(49) 

The maximum velocity is 

v = a{T\ + T 2 ) 

(50) 

and the final position is 

x = + T 2 ){2T\ + T 2 + T^)a 

(51) 

The absolute maxima of acceleration and its derivatives are given by 



\ a ^ max — a M2,k+2( r Tl) , — 0,3 

(52) 


The maxima for the three realizations are compared in table 3. 


Table 3. Absolute maxima for the three realizations of the maneuver M\ 



X 

V 

a 

a 

a 

a( 3 ) 

a o 5) 

300.0 

41.0 

8.7 

4.1 

1.8 

2.4 

( 4 5> ) 3 

300.0 

30.0 

8.2 

3.5 

3.1 

6.2 

(4 5 V 

300.0 

30.0 

7.5 

3.5 

2.7 

7.0 


As noted previously, the position maneuver is smooth but violates the speed limit. The other two 
maneuvers are quite similar and, as shown in figures 20 and 21 , step the aircraft forward 300 feet in 
18 seconds without violating the constraints. It may be noted that, as is to be expected, the satisfaction 
of constraints requires greater activity in the control u. Of course, the increase may be expected to be 
smaller for the optimum solution. 



Figure 20. The three-segment maneuver M 13 = 
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Figure 21. The seven-segment maneuver Mi 7 = 



It is helpful to think of a composite maneuver as a sequence of control cards. For example, the 

three-segment maneuver M\ - (<r ^) 3 may be represented by the sequence C = C 0 C 1 C 2 C 3 of four 
control cards, as follows: 


C 0 = 


Ci = 


c 2 = 


<?3 = 


T = 0 

x = 00000 ** 

T = 8 

x = * 30 0 0 0 14 

T = 2 

x = * 30 0 0 0 14 

T = 8 

x= 300 0 0 0 0 04 


The first card Co specifies the initial condition. The second card C\ commands that the segment is to 
last 8 seconds and that the terminal velocity x = 30 ft/sec (variable 1), higher derivatives up to fourth 
(variable 4 ) are to be zero, while the position is to be ignored. The third card C 2 commands a 2 -second 
coast at 30 ft/sec. The last card commands a steady 300 ft (variable 0). 


This representation becomes particularly helpful when x is multidimensional. Complex, multiaxis 
maneuvers may require many control cards, each with many (e.g., 90) rows of data. In the present 
structure for the flight vehicle management system, these control cards are the interface from the discrete 
levels to the analog levels. Before expanding on this point, consider the multiaxis case. 
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3.4 Coordinate Patching 


To be applicable to aircraft, the preceding discussion must be generalized to include multiple axes 
and rotations. A slight detour into the methodology of dynamic forms is made to considerably simplify 
the discussion that follows. Much more detail is given in reference 3. 

A scalar form of order k of a scalar variable of time is defined by the k-tuple consisting of the 
variable and k of its time derivatives: 

SF k {x ) = x e R (53) 

For example, the system in figure 15 with n = 4 may be represented by SF 5 (x i). Similarly, vector 
and matrix forms are given, respectively, by 

VF k (v ) = (u, v, d( 2 ), . . . , yW oe R n (54) 

and 

MF k (M) = ( M, M, 2 ), . . . , A/W ) , M eR niXn2 (55) 

As described in reference 3, sums, products, and other elementary functions are easily extendible to 
dynamic forms. Then, more-complicated functions may be built out of the elementary set. Thus, for 
example, the possibly time-varying affine map 

v = Mu + z (56) 

may be raised from order zero to order k as follows: 

VF k (v) = MF k (M) * VF k (u) + VF k (z) (57) 

Let a be a three-dimensional, right-handed, orthonormal axis system. The a -coordinates of a 

vector x will be denoted by the column x a . The attitude of another axis system b may be defined 

relative to a by the direction cosine matrix C^ a . The coordinates transform according to x 5 = Ci a x a . 

^ ^ ^ ^ 

Let the angular velocity of b relative to a be denoted by the vector uj ba -, its body coordinates are 
Then the time derivative of the direction cosine matrix, 

^ba ~ S(u}i)at))Cba (58) 

where, as before, the skew-symmetric matrix, given for any z E R? by 

0 23 -22^ 

S(z) = -23 0 21 (59) 

V 22 -Zl 0 
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represents the cross-product operator. The ability to generate direction cosine matrices by a structure 
of the type shown in figure 22 is of interest. It is convenient to represent this structure as the rotation 
form: 


RF k (C ba ) 


^ba > ^ bab > ^bab > ^ 


( 2 ) 

bab ’ ' ' 


u>, 


,(*- 1 ) 
bab 


(60) 


co 


(k-1) 

bab 


co 


babf 


jl _ 


|0) 


bab 



C ba 

J 

C ba 


□ 








Figure 22. Rotational form RF k (C ba ). 


The closely related matrix form is given by 


MF k (C ba ) = 


^ba^ba^ba 


(61) 


The kinematic equation (58) and 

SfWtat) = Cfc,C£ (62) 

where () r denotes the transpose, provide a nonsingular link between the two forms. The subroutines 
for transforming them into each other are denoted as 

MF k (C ba ) = MF^RF[RF k (C ba )\ (63) 

= MF>-RF[MF k (C ba )] 

The generalization of the Coriolis’ derivative to order k is given by 

VF k (v b ) = MF k (C ba ) * VF k (v a ) = * FF^a) (64) 

This subroutine outputs the b-components v b of v and the k time derivatives of v b for given a-components 
v a and its k time derivatives and the attitude C ba and angular velocity and its derivatives. 


The transpose of a matrix form, 

[MF k (M)} T = MF k (M T ) 

Consequently, the inverse of a rotation form is given by 

RF k (C ab ) = MF>-RF{[MF k (C ba )] T } 


(65) 

( 66 ) 


is. 
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It is often necessary to represent rotations in terms of Euler angles in some desired sequence. That 

C b a = E qi ( a l) E q 2 MEq 3 (<* 3 ) (67) 




where the sequence q = (</i, <72> # 3 ) with q 2 7 ^ 71 , 73 , and Ei is an elementary rotation about axis i. 
This representation has a gimbal lock (the Jacobian matrix is singular) at «2 = tt/2 for a nonrepeating 
sequence, and at |c* 2 | = 0 for a repeating sequence. Since derivatives of Euler angles are needed, an 
angle form is defined as follows: 



and the subroutine for the forward and inverse parameterizations (away from the singularities) is defined 
as: 

MF k (C ba ) = MF<AF q \AF k (c ba )} 

AF k (a H ) = MFyAF q [MF k (C ba )\ 

This algorithm provides a bidirectional link between a direction cosine matrix and the corresponding 
angular velocity and its derivatives on the one hand, with the corresponding Euler angles and their time 
derivatives on the other. 

All these easily coded algorithms are very useful in practice. Now consider a maneuver such 
as the transition from VTOL to VATOL, as shown in figure 11. The standard representation by the 
nonrepeating sequence q = (1, 2,3): 

Ctr = £123(7) = £1(71) £2 (72) £3 (73) ( 70 ) 

which is nonsingular for near horizontal attitudes, becomes singular in the vertical attitude. A repeating 
sequence such as p = (3,2,3) 


Ctr = £ 3 .23(0) = £3 (fil) £ 2 (/^2 )£ 3 ) (71) 

is nonsingular near vertical attitudes. To achieve large-angle maneuvers, it is necessary to patch coor- 
dinates smoothly. That is easily done by means of the algorithms defined previously. 

AF^(0) = HP [AF$] = MF^AF{MF<AF [aF*( 7 )]} (72) 


Let the coordinate-designating mode variable be denoted by n a - This mode variable may take on 
many values, but for the VTOL to VATOL maneuver, two are noted: horizontal attitude (HAT) and 
vertical attitude (VAT). Switching from one coordinate system to another is a discrete event initiated by 
a change in fi a . In the horizontal mode fi a = HAT, both sets of coordinates are computed according 
to 


( XFJfr) \ ( AF k ( 7 ) \ 

\AF k (H)) \ H%[AF q {-y)\ ) 


(73) 


That is, 7 is free to move in any preassigned manner while 0 is slaved to the image of 7 . On the other 
hand, in the vertical mode fi a = VAT, 0 is free while 7 is slaved: 


(AF k (-,)\ = ( H}{AF k (f))]\ 
\AF k (0)J ( AFp(0) ) 


( 74 ) 
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The coordinate switch may be commanded by means of control cards. For example, consider a sequence 
of control cards C = . . . CjCj+j . . where 


0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 


0 0 0 0 04 

0 0 0 0 04 

0 0 0 0 04 

0 0 0 0 04 

0 0 0 0 04 

11 0 0 0 14 

Control card Cj commands a segment of 2 seconds duration in the horizontal attitude coordinatization 
mode n a = HAT, in which the 7s are free to satisfy the indicated terminal conditions: zero roll 
(71 = 0), constant pitch rate 72 = 10 deg/sec, and constant yaw rate 73 = 11 deg/sec. The /3s are to be 
slaved to the image of the 7s. The next control card C^ + i commands a segment of 7 seconds duration 
in the vertical attitude coordinatization mode / i a = VAT. The initial condition to be satisfied by the 
(3s are given by the slaved /3s of the previous card Cj. The terminal conditions are constant ft = 0, 
ft = 90, and (3% = 11. The 7s are slaved to the image of the (3s. 

As an example of this type of coordinate patching, consider a large-angle maneuver, as shown in 
figure 23. The maneuver lasting 40 seconds pitches the time-axis system from zero to 90 degrees and 
then back to zero while yawing through the complete 360 degrees. In the vertical position, (t\ ~ —ft), 

t is tilted about £3 from zero to 3 degrees and back to zero. The complete maneuver is synthesized 
with eight control cards. The evolutions of 7 and 7 are shown in the first two panels of the figure. The 
Euler angle forms are set to zero within one degree of their respective singularities. The evolutions of 
(3 and (3 are shown in the bottom two panels. The 3-degree tilt is shown in the top panel of figure 24. 


Ci 


and 


a 


i+1 


/i a = HAT 

71 

72 

73 

Pi 

P2 

ft 

T = 7 

/ia = VAT 

71 

72 

73 

Pi 

ft 

ft 


/ 0 0 
/ 0 10 
/ 0 11 
s 0 C 
s 0 C 
s 0 C 


s 0 
s 0 
s 0 
/ 0 
/ 90 

/ 0 


04 

14 

14 

04 

04 

04 


The corresponding evolution of the direction cosine matrix Cft may be computed as discussed 
previously by means of the algorithm 


f MFyRF{MF^AF[AF^{-/)}} if /i a = HAT 
\ MF>- RF{M F -< AF[AF^(P)] } if /i a = VAT 


The angular velocity, angular acceleration, and their derivatives are shown in in the bottom panels of 

4 — ► 

the figure. The components of uitrt are measures of the rate gyros attached to t . Note that 7, which 

4 — * 

looks like yaw rate, becomes roll rate as t becomes vertical and then again yaw rate at the end of 
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the maneuver. The 3-degree tilt of t\ submaneuver has no effect on the roll channel but activates both 
pitch and yaw. Finally, it may be noted that, as required, both angular velocity and acceleration are 
continuous throughout the maneuver, despite segmentation and coordinate changes. 



t i i i i i i i d 

0 5 10 15 20 25 30 35 40 

time (sec) 


Figure 23. Vertical attitude maneuver: Euler angles. 



0 5 10 15 20 25 30 35 40 

time (sec) 


Figure 24. Vertical attitude maneuver: angular velocity. 


33 





3.5 Transition Dynamics 


In the preceding discussion of the coordinate change 7 <-► (3, there was enough freedom to patch 
continuously across the mode switch. That is not always the case. Two very useful coordinatizations of 
Ctr, namely the runway path axes RPA and wind path axes WPA, are given in terms of the trajectory 

4 — ► i — ► , 

of the origin of t , as follows. Consider an axis system v , defined so that the first unit vector v\ points 
along the velocity v of *v relative to the runway-fixed axes r . In addition, let the total specific force 

f = v-gr . 3 (76) 


be in the (v\,vz) plane, pointing up. That is, 

V\ = v/\v\ 

v 2 = v\ x f/\v\ X /| (77) 

£3 = V\ X V 2 

Then the attitude of *v relative to 7 is given by the direction cosine matrix C vr , whose rows are the 

r-coordinates of v : 

rp 

C V r ~ ( Vlr 1 ^ 2 n ^3r ) (78) 

This zero-order computation can be easily raised to order k by means of the methodology of dynamic 
forms: 

MF k (C vr ) = MF^VR[VF k+l (v r )] (79) 

which may then be transformed into the corresponding rotation form RF k (C vr ). In the runway path 
axis mode (fi a = RPA), MF k {Ct r ) = MF k {C vr ). 


Another frequently used axis system w is defined in terms of the relative wind. The velocity 
relative to the average wind is given by 

v 0, = v — w (80) 

where w is the average wind velocity with respect to the runway. The axis system is given by 


W\ = i7 a /|i? a j 

uf 2 = wi x //|u)i x / 1 (81) 

tog = till x w 2 

so that the first basis vector points in the direction of the relative wind F 1 rather than in the direction 
of v. In the wind path axis coordinatization mode (g, a = WPA), MF k (Ct r ) = MF k (C wr ). Other 
variations on this theme are, of course, possible. The problem is that, while it is possible to switch 
smoothly from a path-defined mode to an Euler angle mode, the change from Euler angles to path may 
be discontinuous because the path is generally not a free variable that can be slaved to an attitude. 
However, the discontinuity may be smoothed by means of transition dynamics, as follows. The idea is 
to factor the attitude into two rotations, one of which is relaxed to the identity. For example, suppose 
that the mode is switched at t — t\ from g, a = HAT to fi a = RPA: 


r _ f Emit), for t < ti 

tr l C vr (t), for t\ < t 


(82) 
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Then, without smoothing, the discontinuity at t\ is 

D (h) = Eiwiti )Cy r {t{) 

On the other hand, if we define 

C ={ ^123(0, for t < t\ 
tr l D(t)C vr (t), for ti < t 


(83) 


(84) 


MF k (C tr ) = 


(85) 


and relax D(t) — > I, then the transition can be made smooth. This zero-order procedure can be easily 
lifted to order k: 

_ f MF k (Ei 2 z), for t <t\ 

~~ 1 MF k (D) ★ MF k {Cvr), for t x < t 

There are many ways to relax MF k (C) — » (I, 0 , One possibility is to express D in terms of 

Euler angles 6 in some sequence q and then drive AFq(6)(ti) — > (0, ...,0) by means of polynomial 
segments. Another way is to apply a zero-seeking regulator structure. For example, suppose that D is 
represented in terms of Euler parameters (e, 77 ) e B? x R : 


and conversely, 


D = I + 2r]S(e) + 2 S 2 (e) 

( 86 ) 

V = b[tr(D) + l]z 
e = ^ axis(D)/rf 

(87) 

axis[S(x)] = x 

( 88 ) 


where for any x G i? , 


The constraint is e 2 + r/~ = 1 , and the representation is singular for |e| = 1 . These zero-order compu- 
tations are easily raised to order k, as follows: 


VF k {e),SF k (r})) = M F >- PF[M F k (D)} 

MF k (D) = MF ~<PF[VF k (e), SF k (rj)] 


(89) 


Associated with the vector part of the parameter form is the three-dimensional string of k integrators, 
each 


The initial condition is given by 


Choose a stable regulator law 


constructed as follows 


MF k (D) 


e ( fc+1 >=« 

(90) 

,..,eW(t 1 )]=VF k (e)(t l ) 

(91) 

u = k e, . . . , 

(92) 

will decay to zero. The evolution of Ctr f° r t\ 

< t is then 

- [1 - VF k {e ) 2 ]2 

= MF<PF{VF k (e),SF k {7])} 
= MF k (D) * MF k (C vr ) 

(93) 
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By construction there is no discontinuity in Ctr or k of its time derivatives, and Ctr coincides with C vr 
after the transient has died down. 

Several useful ways have been developed to coordinatize large-angle attitudes. The generation of 
three-dimensional reference trajectories is now considered. 

3.6 Output Mode 

The coordinate mode, /i a , discussed in the preceding section, affects the meaning of the system 
state and control coordinates. In this section the output mode, n y , which determines the meaning of the 
system output (that is, which variables are to be tracked) is considered. Some of the many possibilities 
are the Cartesian, cylindrical, or spherical coordinates of position, or airspeed and path, or, as in the 
360-degree roll maneuver shown in figure 14, speed and lateral and normal forces. Many more output 
maps are of practical interest. Specifying and switching of output maps may be handled in much the 
same way as the Euler angles. Thus, for example, consider a 90-degree horizontal turn maneuver with 
100 ft/sec initial and final speeds. The maneuver may be realized by means of the Cartesian coordinates, 
as follows: 


M X y’. 90-degree Cartesian maneuver 


T = 0, (ia = RPA. fiy = CAR 

VF^Tr) f 
VF 2 (v r ) f 
VF 3 (r r ) f 

T = 16.5, Ha = RPA, Hy = CAR 

VF^Tr) f 
VF 2 (r r ) f 
VF s (r r ) f 


0 100 0 0 0 
0 0 0 0 0 
0 0 0 0 0 


1000 0 0 0 0 04 

1000 100 0 0 0 04 

0 0 0 0 0 04 


The turn is accomplished with one segment in 16.5 sec. The path is shown in figure 25. The attitude 

of t is slaved to the runway path axis system (/x a = RPA). The corresponding Euler angles a 
in the standard sequence q = (1,2,3) may be extracted from Q r . Figure 26 shows the Euler angles, 

i ^ ^ y 

angular velocity, and angular acceleration of the t system relative to the runway r . Figure 27 shows the 
evolution of the ground speed v (ft/sec, upper panel), acceleration v (ft/sec 2 , second panel), f-coordinates 
of force f[ (third panel), and curvature k (in inverse thousand feet), where 

v = 
v = 

ft = 

K = 


t[ r v r 


Ctrifr + 963 )/ 9 
\\d/dt(^)\\/v 


(94) 
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time (sec) 


Figure 27. Cartesian maneuver: v (ft/sec), v (ft/sec 2 ), f[ (g’s), k (1/kft). 


and where g is the acceleration of gravity. It may be noted that the ground speed of this maneuver is 
not constant; there is a small acceleration into the turn and a small deceleration out of the turn (second 
panel). In some cases a constant speed maneuver is preferable. That may be accomplished by means 
of the cylindrical coordinates (v,tp,h) of velocity: 


v r = 


v cos ip') 
v sin ip 
~ h 


(95) 


where h — —r r 3 is the altitude. The Cartesian coordinates of position and the arc length are, of course, 


r r (f) = r r ( 0 ) + / v r (a)da 
J 0 

s(f) = s( 0 ) + J v(a)da 


(96) 


Often it is easier to generate the maneuver in two steps. First the path is “drawn” with some convenient 
speed, and the path parameter is changed from time to arc length. Then the path is traversed at some 
other speed. This may be accomplished as follows: 


First draw the path, as a function of time 

(s,^) = [si(f),^i(f)], t £ [0, T] (97) 

where sj > 0. Next consider s as the independent parameter: 

(t,^) = [a(a), /?(«)], *€[*i(0),«i(T)] (98) 

where si[a(s)] = s, and p(s) = '0i[a(s)]. Finally, traverse it with a different schedule 

s = S 2 {t), t £ [t 0 ,tf] (99) 
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namely 


(s,ip) = [S2(t) t ih(t)] = {-*2(0* >®[«2(*)]}» t € [t 0 ,t f ] (100) 

where to = a[si(0)] and tf = a[si(T)]. The problem is to raise this zero-order computation to order 
m. That is, the problem is to construct 

SF i W)=V>l{a[SF t (s)]} ( lOt) 


Suppose that for given sj = ao the corresponding t\ and SF k (s\) are known. Then the evolution 
of sj near the base point (ao,ti) is given by the Taylor series in r, 

si(*l + r) = + sir + . . . -1- sf\ k /k\ (102) 

This is a polynomial of degree k in r with constant coefficients 

a ~ (si,si,...,s[ k ^/k\) (103) 

This zero-order polynomial may be raised easily to order k: 

SF k (si) = poly[a,SF k (r)] (104) 

The input is SF k (T\) = (r, 1,0, . ..,0) and the output is SF k (s\). The present objective is to invert 
the process: given some other evolution SF k (s2), find the corresponding SF k (r2). Having found that, 
the new evolution of ip is given by the polynomial 

SF k (ip 2 ) = poly[b,SF k (T2)] (105) 

with constant coefficients 

The inverse SF k (r2) may obtained in the standard way. (See ref. 3.) 

SF k (r) = (0, . . . , 0) 

9 = (s) _1 
do i = ljk 

w = {poly [a, SF k (T)]} k 
r( k ) = g(s ( k ) — w) 
end do 

In the special but useful case when both the initial drawing and final traverse are done at constant but 
possibly different speeds, si = v\ and «2 = v 2 = respectively, simple scaling suffices: 

Sf*(V*) = (0 i.w/’i W*) C07) 

For example, consider the following abbreviated control card schedule. 
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M v ^\ 90-degree cylindrical maneuver 


0,/J y = CYL 

SF 4 (s) 

f 

0 

100 

0 

0 

0 


SF 4 {^) 

f 

0 

0 

0 

0 

0 


SF 4 (h) 

f 

100 

0 

0 

0 

0 


5, \x y = CYL 

SF 4 (s) 

f 

0 

100 

0 

0 

0 

14 

SF 4 (ip) 

f 

0 

6 

0 

0 

0 

14 

SF 4 (h) 

f 

100 

0 

0 

0 

0 

04 

10, fJ,y = CYL 

SF 4 (s) 

f 

0 

100 

0 

0 

0 

14 

SF 4 {iP) 

f 

0 

6 

0 

0 

0 

14 

SF 4 (h) 

f 

100 

0 

0 

0 

0 

04 

5, n v = CYL 

SF 4 (s) 

f 

0 

100 

0 

0 

0 

14 

SF 4 (xp) 

f 

90 

0 

0 

0 

0 

14 

SF 4 (h) 

f 

100 

0 

0 

0 

0 

04 


The speed is held constant throughout the maneuver at v = 100 ft/sec. At the end of the third segment, 
the displacement r r (20)— r r (0) = (1208, 1208, 0) T . A displacement of (1000, 1000, 0) may be obtained 
by scaling time by a factor of 1000/1208. The resulting path is shown solid in figure 28, and it may 
be compared with the Cartesian path (dotted). The behaviors of the attitude, angular velocity, and 



Figure 28. Cylindrical maneuver: path. 

acceleration are shown in figure 29. Figure 30 shows that the speed (top panel) is constant, as expected. 
The curvature is smaller for this case than for the Cartesian maneuver shown in figure 27. 
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Often ground speed is of lesser interest than the airspeed. The air velocity is given by 

v a = v — w = st — w (108) 

where t is the path tangent. Therefore, the ground speed needed to maintain a desired airspeed v a = a, 
which may be a function of time, is given by 

s = t ■ w + \J (t • tu) 2 + a 2 — iu 2 (109) 

This equation can be used iteratively to raise the order to SF k (s). Often it is desirable to point the t 
axes into the wind. This maneuver may be done by choosing t\ in the direction of the air velocity v 0, 
rather than, as before, in the direction of the ground velocity v. 

Figure 31 shows the 90-degree turn in which the airspeed is held constant at v a = 100 ft/sec in the 
presence of a constant wind w T = (—20, 20, 0) ft/sec. The maneuver is executed in the wind path axis 
mode fj, a = WPA. It may be noted that, while the ground speed v (solid line in panel 1) varies, the 
airspeed v a (dotted) is constant at the set value of 100 ft/sec. The second panel shows f[. It may be 
noted that only the normal component (dashed) is active. The Euler angles at r degrees are shown 
in the third panel. Note that, since the axis is pointing into = v r — w r , the yaw angle atrz (dashed) 
starts the turn with —11.5 degrees and finishes the turn with 78.6 degrees. Furthermore, the roll angle 
(solid) is not constant in the constant radius segment. The bottom panel shows the angular acceleration 
u>trt degrees/sec 2 . Note the smaller roll acceleration at the beginning of the maneuver and the larger 
value at the end in comparison with the no-wind case shown in figure 29. Finally, the duration of the 
maneuver is a little longer, namely 17.35 seconds. 



time (sec) 

Figure 31. Fixed air speed cylindrical maneuver {v,v a ),f[, at r ,d>trt- 
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Other possibilities are of practical interest. For example, spherical coordinates of velocity 
= SPH) 

v r = ( v cos if ) cos 7, v sin ip cos 7, — v sin 7 ) T (11 0) 

are frequently used for guiding an aircraft to the airport. 

The methodology developed in the present chapter allows one to generate segments in any one 
of these many useful representations, and to patch the segments into a smooth, executable maneuver. 
Thus, the Guidance Output Trajectory Generator, which is the first block in figure 5 , is realized. The 
input is a sequence of control cards. The output is the commanded definition and evolution of the 
system output y c and the evolution of the operating mode So far, the operating mode had two 
coordinates, namely y a , which defines the coordinatization of the system state space, and y y , which 
defines the state-to-output map. In the next chapter, in which the Guidance State Trajectory Generator 
is considered, the third coordinate, the control mode /r c , which specifies how the maneuver is to be 
executed, is introduced, and an approach for constructing the guidance state and control ( x 9 , u 9 ), which 
produce the desired output y 9 ~ y c , is given. 
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4 DYNAMIC INVERSE 


In the previous section, a procedure for converting control points into the desired motion of the 
output is developed. (The underline denotes a function of time.) In this section, an inversion 
procedure that converts y ^ into the corresponding motions of the state and control (x 9 ,u 9 ) is developed. 
It will be useful to clearly distinguish between two types of inverses, namely the static inverse and the 
dynamic inverse. 


For the static inverse problem, the following terminology is employed. A direct algorithm / is 
given: 

u e U c R m , p e P c R q , y g Y c R m 

f : UxP-^Y (111) 

y = f{u,p ) 


where u, p, and y are, respectively, the control, parameter, and output. The problem is to construct the 
relative inverse algorithm /, 


/ : Y x P-^U 

u c = f(y<bP ) 


( 112 ) 


so that, for any desired output y^ G Y and parameter p G P, the output error 


is sufficiently close to zero. 


e y = Vd~ f(u°,p) 


(113) 


In contrast, the dynamic inverse problem may be summarized as follows: Suppose that a state 
equation and an output map are given: 

t e It c R, X € A c R n , U e U c R m , y g Y c R m , p e P c R 9 
f : X xU X P X 

h : XxP^Y (114) 

x = f(x,u,p ) 
y = h{x,p) 

Here, x, u, p, and y are, respectively, the state, control, parameter, and output; / is the flow, and h is the 
output map. In this case, the inversion problem is to construct an algorithm T, which transforms any 
given, suitably restricted, motion of the output y^ and parameter p, respectively, into the corresponding 
motion of the state x 9 and control u 9 , which satisfy both the state equation and the output map. That 
is, the state error equation and output equation error, respectively, 

e x (t) = ± 9 (t) - f[x 9 (t),u 9 (t),p(t)] 
ey(t) = Vd(t) ~ h[x 9 (t),p{t)} 

must be sufficiently close to zero on the interval of interest /<. 


In the case of linear systems, the dynamic inverse is given by the particular solution. For example, 
consider a linear, time invariant, scalar, n th 

-order system, 
x G R n , u G /?, y G R 

x = Ax + bu (116) 

y = cx 
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As is well known, if the system is controllable and observable, then there exists a transformation of 
state and control coordinates 


2 = Tx 
v — wu + rx 


(117) 


with nonsingular T and w, which turns the state equation into a simple string of integrators (Brunowsky 
form with the Kronecker index k = n) 


z k = z k+ 1 » 1 < k < n - 1 
z n = v 

For the system considered here, the rows of T and the row vector r are given by 

t( — fi-A* -1 , 2 <i<n 
r = tiA 71 

while the first row of T is found from the n equations 

t\A l ~ l b = 0, l<i<n— 1 
t\A n ~^b = w 

Then, for cT~ l = (ho, , h n -\), the output is given by 

n_1 tu\ 

y = E h k z[ k) 

k = 0 

But for a given y^, that is a driven differential equation 

h k z \ k) = Vdit) 

n — 1 


(118) 


(119) 


(120) 


( 121 ) 


( 122 ) 


This equation is referred to as the driven zero dynamics of system (116). Its general solution is given 
by the sum of the particular solution z\ p and the homogeneous solution z^. Rather than the solution 
of the usual initial value problem, the “steady-state” (bounded-input-bounded-output) solution of this 
equation, which is just z\ p , is desired. Thus, if the desired output is 


y d (t) = s\n(ujt), teR (123) 

then, of course, the particular solution is given by 

zi p (t) = A sin(cut + 0), t € R (124) 

where A and 0 are the gain and phase of the transfer function, 

AeP^ = g(ju>) = h(ju)~^ (125) 


The roots of h( A) = 0 are assumed to be far enough from the jtu-axis: |i?e(A)| > a for some given 

a > 0 . 
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The particular solution for the more general case of the forcing term, y^, may be solved by means 
of a convolution. Let g(t) denote the Fourier inverse transform of the transfer function, g(juj). Then 

/ oo 

9(t - r)y d (T)dr (126) 

-00 

In practice, the integral may be approximated by 

zip(t) = j_ a g{t - r)y d {T)dr (127) 

where a is only several times the longest time constant l/a of g. 

The corresponding state and control [z p (t), v p (t)] are given by the appropriate derivatives of z\ p (t), 
and the solution to the inverse problem is given by 

x 9 {t) = T~ l z p {t) 

u 9 (t) = w~ l [v p (t) - rT~ l z p {t)\ 

Note that, since the homogeneous solution is chosen to be zero, it does not matter whether the differential 
equation (122) is stable or not, that is, whether the zeros dynamics are or are not minimum phase. On 
the other hand, the computation may be noncausal with respect to y d , requiring a “look-ahead” at t of 
a units of time: 

Vd(t + r), —a < t < a (129) 

The computation of the particular solution for the multiaxis case ( u , y € R m ) is the same if the 
Brunowsky form is rectangular, that is, if all m strings of integrals are of the same length n = n/m. 
Otherwise, the computation is similar but a bit more complicated, requiring some additional bookkeeping. 


In many cases, it is not necessary to obtain a particular solution that satisfies equation (122) 
because often the zero dynamics are caused by small, essentially parasitic, effects so that the zeros of 
equation (122) are far away from the spectrum of the forcing term y d . In this report, the cases where 
the ho term is nonsingular and dominant are of interest, while the rest of the terms represent a small, 
parasitic perturbation. This distinction is made explicit by means of the parameter e: 

h{D,e) = eihn^D 71 - 1 + . . . + hxD 1 ) + ho (130) 


For the case of e = 0, the particular solution is just 

zip(t) = hy l y d (t) (131) 

Then, in view of equation (118), and because ho is constant, ( z p , v p ) may be obtained by differentiation 
of the desired output y d : 


= h o 1 y d k l \t), l<k<n-l 

= WW 


The corresponding evolution of the natural 
equations (117); that is. 


Xp(t) 

v?p{t) 


state and control coordinates (x°, u°) is obtained by inverting 


= T~ l z p (t) 

= w~ l [v p (t) - ra;jj(0] 


(133) 
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This solution is referred to as the pure feedback approximation. Suppose that 

x 9 (t) = x°(t) 
u 9 (t ) = u°(t ) 

Then the state equation and output errors are given by the linear version of (115) 


( 134 ) 


e x (t) = x 9 (t) - [Ax 9 (t) + bu 9 (t)\ 
e y(t) = Vdi 1 ) - cx9 (t) 

In the absence of zero dynamics, that is when = 0 for k > 0, the pure feedback approximation is 
exact in that it satisfies both the state equation ( e x = 0) and the output ( e y = 0). Otherwise, there will 
be errors in either one or the other or both equations. Often in practice these errors are completely 
negligible. However, when the zeros are closer in and the desired maneuver is relatively fast, the errors 
due to the neglected zero dynamics may become significant. In that case, the zero dynamics must be 
accounted for. 


This approach may be extended to nonlinear, time-varying, multiaxis systems as described in 
references 7-10. As in the linear case, it is useful to introduce e into the state equation and the output 
map 


x = f(x, u, e) 
y = h(x,e) 


(136) 


so that the zero dynamics, if any, may be suppressed from the model by setting e = 0. The pure 
feedback solution (x°, u°) obtained for the model with e = 0 may be taken as the first approximation 
of the inverse. The errors due to the neglected dynamics are given by 


4(0 = x°(t) - f[x°{t),u°{t),l\ 

4(0 = Vd( t)-h[x°{t),l\ 

An iterative (Picard) procedure for reducing these errors is described in reference 9 and applied in 
reference 10. The procedure is based on the convolution (127). In the present report the computation 
of the pure feedback approximation, and controlling the resulting errors, is of interest. 


A very easy and practical way to satisfy the state equation so that e x = 0 follows: Consider the 
servo structure shown in figure 32. The left block represents the approximate inverse T° such as the 
pure feedback algorithm in which the zero dynamics are suppressed (c = 0), or the pure feedback 
approximation corrected by only one Picard step. The right block represents the simulation of the exact 
model (e = 1). The input to this block is the control u 9 \ the output is the state x 9 . So, by construction, 
e x = 0. The tracking error ey is regulated by the middle block. That is, a regulator law 


u gr = k(x°, x 9 ) 


(138) 


is designed and the loop is closed: 


u 9 = u° + u 9r (139) 

Then it may be expected that, for a reasonable regulator, sufficiently weak zero dynamics, and suitably 
restricted class of inputs y^, the tracking error will be small. 
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guidance regulator 


Figure 32. Servo structure satisfying the state equation. 


Next, the key problem of computing the pure feedback approximation is considered. 

4.1 Dynamic Inverse of Pure Feedback Systems 

In the preceding section, the dynamic inverse of a linear system (116) was considered. The key steps 
in the construction of the pure feedback approximation were the static inversion, as in equation (131), 
and the differentiation of that inverse, as in equation (132). The computations may be generalized to 
nonlinear cases as follows: Consider again a function 

x € X cR m , pePcR?, y<EY cR m 

f : XxP^Y (140) 

y = f(x,p) 

and suppose that this (zero-order) algorithm / has been raised to arbitrary order k, which is a key step 
in what follows. The corresponding algorithm is denoted as follows: 

VF k (y) = f[VF k (x),VF k (p)\ (141) 

That is, equation (141) is an algorithm whose input consists of x and k of its time derivatives, namely 
the vector form 

VF k (x) = (142) 

and similarly for the parameter p. The output of the algorithm is y and k of its derivatives 

VF k (y) = [y,y,...,yW] (143) 

Inversion of this algorithm is desired, so that for given evolution of the output VF k (y) and parameter 
VF k (p) the new algorithm produces the required input VF k (x). 

Now, suppose that for / in equation (140), the Jacobian matrix f x is nonsingular in the region of 
interest, and that the inverse algorithm / 

/ : Y x P —> X 

x c = fiVdiP) (144 ) 

e v = yd-f( xC >p) = 0 
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is available. This zero-order inverse may be raised to order k, as follows. For k = 0, the zero-order 
inverse (144) computes the base point: 

x = f(y,P ) (145) 

At that base point, the first derivative of y is given by 

V = fxi + f p p= {/[(*, 0), (p,p)]}i + f x x (146) 

where (x, 0) denotes the vector form of order one in which the first derivative is zero. The subscript 1 
denotes the first derivative of the variable, namely, 

{VF k (y)} l = y (147) 

and it is evaluated by the algorithm given in equation (141), the vector form of f(x,p). Similarly, for 
the second derivative 


j) = fxi + f x x + fpp + fpp = {/[(x,x,0), (p,p,p)]} 2 + fx x (148) 


and so on for higher derivatives. Hence, the following algorithm for passing derivatives through an 
inverse: 


x = f{y,p) 

VF k (x) = (x, 0, . . . ,0) 


do j = 1 , k 

= fx 1 



end do 


{{[VFi(x),VFl( V )\) J \ 


(149) 


It should be noted that, remarkably, only the base point x and the Jacobian f x at the base point are 
needed to generate the derivatives of x from the derivatives of y and p. 


Next, consider the construction of a dynamic inverse as in equation (114) for the special case of 
pure feedback systems. Later in the present report it will be shown that, for the purposes of guidance, 
the following model may be used to approximate many aircraft types. (See figure 36.) In this model, 
the structure of the multiaxis state, control, and parameter is given by 

Xi,X2 € R ni , X3,X4 ,u € R n2 , p G R np (150) 

with n 2 > ni . The state equation is of the following form with sufficiently smooth flows /j. 


X\ = /i(a?i, 

x 2i 

0 , 

0 , 

0 , 

p) 

x 2 = f2( x l, 

x 2i 

^3. 

0 , 

0 , 

p) 

x 3 = M x l> 

x 2, 

•^3? ^4? 

0 , 

p) 

x 4 = h(xi, 

X2, 

^3? 'E'A') 

u, 

p) 


The structure of the multiaxis output is given by 

yi e R n \ 7/2 e R n2 ~ ni 


with the sufficiently smooth output map 

0, 0, p) \ 

0, o, p ) ) 


(yi \ = ( h i( x i i °> °> 

\7/2/ V ^2(^1 > 


(151) 


(152) 


(153) 
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The block diagram of the model is shown in the upper part of figure 33. The heavy lines are no- 
dimensional. Note that, disregarding the parameter p, there are no lines bypassing integrators in the 
forward direction. Hence there are no zeros. There is a new (canonical) set of coordinates of state and 
control related to the natural state and control ( x , u) by the nonsingular transformation T, in which the 
system takes the particularly simple (Brunowsky) form shown in the bottom of the figure. There is a 
bundle of n\ strings of integrators, each four integrators long, and another n 2 — ni dimensional bundle 



Figure 33. Pure feedback model of aircraft. 


of strings, each of length two. The inverse transformation T converts the Brunowsky form back into 
the natural representation. The construction of the dynamic inverse is equivalent to the construction 
of T : the evolution of the output \VF^(y\), VF^{y 2 )\ and parameter VF^(p) is transformed into the 
corresponding natural state and control (x,u). Both the inverse f and the direct T transformations are 
constructed as follows: 

For the construction of T, assume that VF k {y\) and V F k ~ 2 (y 2 ) are given and that all the functions 
and their inverses have been raised to sufficient order, as discussed previously. Then from the inverse 
h\ of hi, the derivatives of x\ are: 

VF k {xi) = hi[VF k (yi), 0,0, 0, 0, (154) 

and consequently k — 1 derivatives of xj, namely VF k ~^(x\). Then, inverting the first line of the state 
equation (151) gives the derivatives of X 2 '. 

VF k -\x 2 )=h [VF k -\ Xl ),VF k - l (xi),0,0,0,VF k - l (p)] (155) 
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and VF _2 (£ 2 ). Next, since there are ni equations from / 2 and n 2 - equations from h 2 , the 
determination of the n 2 elements of x 3 requires the simultaneous inverse of (/ 2 , /i 2 ) with respect to x 3 
to produce the derivatives of 0:3: 


VF k -\x 3) = UM [VF k - 2 (x l ),VF k - 2 (x 2 ),VF k - 2 (x 2 )A 0 ,VF k - 2 (p) 
and VF k ~^(x 3). Similarly for x 4 and u: 


(156) 


VF k ~\x A ) = hlVF^M, VF k ~\x 2 ), VF k ~ 2 (x 3 ), VF k ~ 3 (x 3 ),0, KF*" 3 (p)] 

VF k -*(x 2 ), VF k ~\. x 2 ), VF k ~ 4 (x 3 ), VF k ~ 4 (x 4 ), VF k ~ 4 (p )] 


At this point, the function 


(157) 


f : (R ni ) k+1 x (R 712 - 7 ^)^ 1 x (R n P ) k+1 -> (i? ni )^ x (i ? n2 )^ -1 x (R n2 ) k ~ 2 x (R n2 ) k ~^ x (R ni } k ~ 4 

(158) 

which maps the desired output and some of its derivatives into the corresponding state and control, and 
some of their derivatives has been constructed: 


T 


VF K ( yi ),VF k - 2 (y 2 ),VF K (p ) 


VF k (x 1), VF k ~ 1 (x 2 ), VF k ~ 2 (x 3 ), VF k ~ 3 (x 4 ), VF k ~ 4 (u) 

(159) 


The multiaxis dynamic inverse algorithm T with k = 4 that is desired is contained in f . The state and 
control (x, u ) are given by the zero-order entries of the output of T, namely. 


(x,u) = T{VF i (y 0, VFHvi), VF 4 (p)] = {t [VFV), VF 2 (y 2 ), KF 4 (p)] } Q (160) 


It is also very useful to be able to go the other way, that is, to have the direct dynamic map T which 
transforms the state and control (x, u) into the Brunowsky canonical coordinates VF 4 (y \ ) and V F 2 (y 2 )- 
This algorithm may be constructed as follows: From the given state and control [VF°(x), l^F°(u)], 
parameter VF 4 {p), and the state equation is obtained 


then 

then 

and finally. 


VF\ Xi ) = fi[VF 0 (x l ),VF 0 (x 2 ),VF\x 3 ),VF 0 (x A ) > VF 0 (u),VF°(p)}, i = 1,4 

VF 2 ( Xi ) = f i [VF 1 (x 1 ),VF 1 (x 2 ),VF l (x 3 ),VF l (x 4 ),0,VF l (p)}, i = 1,3 
VF 3 (x i ) = f l [VF 2 (x l ),VF 2 (x2),VF 2 (x 3 ),0,0,VF 2 (p)l i= 1,2 
VF\ Xl ) = f l [VF 3 (x l ),VF 3 (x 2 ), 0,0,0, VF 3 (p)} 


The outputs are given by 


VF 4 (y i) = h x 
V F 2 (y 2 ) = h 2 


V F 4 (x 4 ), 0, 0, 0, 0, V F 4 (p)] 

VF 2 (x i), VF 2 (x 2 ), VF 2 (x 3 ), 0,0, VF 2 {p)\ 


(161) 

(162) 

(163) 

(164) 

(165) 
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Thus, the function 

T : (fl ni ) 2 X (i? n 2) 3 X {R n P) 4 -> {R n l) 4 X (f? n 2-«l)2 x (/J«l)4 x (fl"l)3 x (#"2)2 x R n 2 (166) 

which transforms the natural state and control (x,u) and the evolution of the parameter VF 4 (p ) into 
the output and state rate and some of their derivatives, has been constructed. 

T(x,u, VF 4 (p ) ) = ( VF 4 (yi), VF 2 (y 2 ), VF z {x i), VF 2 (x 2 ), VF l {x 3 ),VF°(x 4 ) ) (167) 

The transformation T from the natural state and control (x, u ) to the canonical control is just the 
projection P\ 2 of T on the first two slots of its output; 

VF 4 (yi),VF 2 (y 2 )] = T [x,u,VF 4 (p)\ = P 1j2 T [x, u, VF 4 (p)] (168) 

It should be noted that the time derivatives of the outputs y\ and y 2 are obtained without differentiating 
either the state x or the control u. Both the forward map T and the inverse T differentiate known 
functions; neither differentiates signals. 

The preceding discussion may be summarized by means of the block diagram in figure 34, as 
follows. The system x = f(x,u,p), which appears multiaxis, nonlinear, and time-dependent in its 
natural coordinates (x,u), becomes a simple collection of decoupled strings of integrators when viewed 
through the coordinate transformations T and T. 


VF 4 (y 1 ) 

VF 2 (y 2 ) 


VF 4 (p) 

VF 4 (yi) 

VF 2 (y 2 ) 



Figure 34. Linearization of pure feedback systems. 


In the following discussion an aircraft model that will be considered basic for the purposes of 
guidance is developed. Next the notion of aircraft control modes is introduced. Then a procedure for 
computing the pure feedback solution ( x ° , u°) for any of the many possible control modes is presented. 
Finally, a general structure of the regulator is developed. The discussion is illustrated by applications 
to the Harrier aircraft. It should be noted that the procedure is general and applies to many types of 
aircraft and helicopters. 
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4.2 Aircraft Model for Guidance Purposes 

The model developed in the second chapter is summarized and extended as follows. As before, let 

the runway-fixed axes and body-fixed axes be denoted by r and b , respectively. Let r r and v r denote 
the runway coordinates of the aircraft position and velocity relative to the runway, respectively. Let the 
attitude of the aircraft relative to the runway be denoted by the direction cosine matrix C hr , and let the 

body coordinates of the angular velocity of b relative to r* be denoted by u brb . The structure of the 
rigid-body state is chosen to be 


Xw € Xu c ( R 3 ) 2 X [50(3) X 5 3 ] 
Xu) ( T r i v r i C br , UJ b j~ b ) 


069) 


The structure of the control is 

u € U C (R 3 ) 3 

U = ( U m ,U p ,Uf ) 

where u m , u p , and Uf are the moment control, power control, and configuration control, 
there are the parameters with the following structure and interpretation. 


(170) 
In addition, 


p € P C R 3 x Rx R 

p = (t v r ,W,0 T ) (17,) 

where w r represents the wind in runway coordinates, W is the aircraft weight, and 0 T is the temperature 
of the air. The runway and body coordinates of the relative air velocity are, respectively, 


— 

_ 


V r — w r 

C br v r 


(172) 


The input to the force and moment algorithm has the following structure and interpretation 


x b e X b C U x Xu x P 

x b = ( u, Xu, p ) = ( u m , Up, uj, r r , v r , Cfrj., u) brb , w r , W, Bj ' ) 

The output will be expanded to include several other variables of interest: 


(173) 


Vb € Y b = (R 3 f 

Vb = ( m b ,m r ,f b ,a r ,v % ) 


(174) 


In the interpretation, the entries are body and runway coordinates of the generated specific moment, 
body coordinates of the specific force, runway coordinates of acceleration, and the body coordinates of 
the air velocity. Note that the coordinates of y b are not independent. Ultimately, an independent set 
will be associated with each operation mode. 


The algorithm specific to a particular aircraft is the function 


F b : X b -+Y b 
Vb = F b [x b ) 


(175) 
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( 176 ) 


On occasion the three-dimensional (3-D) components of the algorithm F b , 

F b = ( Fl nb , Ff ) nr , F[\ Ff, Ff ) 

will be needed. Next a pair of integrators are added in front of power and configuration controls (up,uj) 
in order to enforce continuity in these controls and their rates. The result is the aircraft model shown 
in figure 35. 



Figure 35. Block diagram of the basic aircraft model. 


The model is quite an accurate and general representation for the purposes of designing aircraft 
guidance systems. The representation of attitude by means of the direction cosine matrix C br is not 
limited to small angles. On the other hand, the presence of the direction cosine matrix makes the state 
space non-Euclidean. Furthermore, for certain aircraft maneuvers, and in certain patches of the state 
space, it is necessary to represent the aircraft attitude in terms of Euler angles in some given sequence. 
Now such state space coordinate changes are considered. 

4.3 Intermediate Axis System 

( > 

The usual procedure is followed; an intermediate axis system, denoted here by f , is introduced, 
and the aircraft attitude, C br , is factored into two rotations, 

C br =C bt C t r (177) 

where Ctr may be a large-angle rotation, but C bt is small enough so that the following Euler expansion 
may be always assumed to be nonsingular: 

Cbt = E 2{^bt\)^{ a bt2) E l(o!bt3) O’ 78 ) 


54 



Several examples of parameterizations of Ctr have been given in section 3. There the mode variable 
fi a was introduced and the horizontal attitude mode HAT, vertical attitude mode VAT, runway path 
axes mode RPA, and wind path axes mode W PA were discussed. 

With respect to the small-angle rotation C bt , the angles a^, and a^ 3 are pitch, yaw, and 

4 — ► i — > i — ► < — ► 

roll angles of the b - axes relative to the t -axes. When t and b are selected to be the wind axes and 
body axes, respectively, these angles become the angle of attack, negative of the sideslip angle, and the 

4— ► W 

roll angle of b . 


Whatever the interpretation of the two factors of C br , the force and moment generators are needed 
in terms of the Euler angles a bt of the small rotation C bt instead of the composite rotation C br . 


The Euler angles a bt and direction cosine matrices C br and Ctr as well as their derivatives to order 
k are easily related to each other by means of the dynamic forms. The angle form 

= (<%,<% (l79) 

and the rotation form 

RF t k r = ( Ctr , u, tTt 4ri l) ) < 180 > 

may be converted into RF^ r by first converting the angle form into the corresponding rotation form by 
means of the algorithm RF -< AF and then multiplying by RF^ r : 


RF&. = RF -! AF,(AF&) . Rff r 


(181) 


Conversely, suppose that RFfL and RF t k r are given. Then the corresponding angle form is given by 


AF& = RF >- AF„ * RF k t ] 

where RF^ t is the rotation form for C r t — Cj r , the transpose of Ctr- 
Consider now a new representation of aircraft rigid-body state. 

x Q € X a C (7? 3 ) 4 


X a 


( 1 ) 

r r ,v r ,a bt ,al t 


(182) 


(183) 


which is obtained from x w by replacing by (a bt ,a bt ). The two sets of variables are related 

by the special case of the algorithm described previously: 

U(,4« l> ) = RF>AF q [RF^,RFt t \ 

( Cj lT . ^7,,./, ) = RF < AFg(AF^ t ) * RF^ r 

The main advantage of the new state space X a is that it is (locally) unconstrained, whereas the old state 
space X u has the orthogonality constraint C br C^ r = 7. 
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It is also desirable to expand the force and moment generator output as follows. 


yt e Y t c (f? 3 ) 7 

Vt = ( m b ,ma,m r ,ft,fb, a r,v % ) 


where m a is the second time derivative of a b t, and ft is the specific force in the t coordinate system, 
namely, 

m a = a = { RF AF 23 i [RF^ * flF,?j]} 2 (|86) 

ft = clh 

Note that, whereas m b and m r represent the specific moment vector, m a does not. Also, the subscript 
231 represents the Euler sequence in equation (178). The new model of the force and moment model 
algorithm will be denoted by Ft: 


x t e X t C U x X a x [50(3) x ( R 3 ) 2 ] x P 

x t = ( u,x a ,RF? r ,p ) 

F t : Xt^Yt 

yt = F t {x t ) 


(187) 


The term in square brackets is the domain of the second-order dynamic form RF^ r of the intermediate 
axis system. The (3-D) components of Ft will be grouped and denoted as follows. 


F t = ( F[ nb , F™ a , F/ nr , F{\ F t fb , F t ar , F t vb ) 


(188) 


The rigid-body state equation in the new coordinates becomes 


fy = v r 

v r = F ( ar (u,x a , RF ? r , p) 
(1) 

a bt = % t 

= Fr(u,x a ,RF? r , V ) 


(189) 


The system is shown in block diagram form in figure 36. The computations of the body direction cosine 
matrix C^, body angular velocity, and acceleration RF^ r are also shown in the diagram. 


The following list summarizes the interpretation of the variables: 


u 

— (%! 

U P> 

Uf) 






Xu 

= (r r , 

v r , 

C b r , 

^brb) 





X a 

= (r r , 

Vr , 

a bti 

q.W') 
a bt ) 





X b 

= ( u m, 

Up , 

Uf, 

r r , 

V r , 

&br > 

UJbrbi 

p) 

Xt 

— ( u m, 

Up ^ 

Uf, 

r r, 

v r , 

a bt> 

a (1) 
a bt > 

Ctr, LjJ trt 

yb 

= ( m b , 

m T) 

fb, 

Clr, 

4) 




yt 

= (mb, 

m a , 

m r 

ft, 

f b , 

Ur ? 

4) 



&trti P ) 


(190) 
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Figure 36. Alpha guidance model of the aircraft. 


In a simulation, the initial conditions on the integrators in figure 36 and an evolution of all the inputs 
are given, and the simulation output is the evolution of all the states and various functions of the states. 
In a typical guidance problem, the evolution of the aircraft position r r is given, and the problem is to 
fill in all the rest of the variables in the diagram. There are many ways to do that. Just which of the 
solutions is actually wanted and what additional information is to be provided is defined by the control 
mode. This mode is discussed next. 


4.4 Control Modes 

Now, the objective in this report is to construct pure feedback approximations of the inverses of 
the aircraft equations of motion. Consider, again, the model shown in figure 35. Suppose that at each 
instant of time t the desired evolution of aircraft position VF*(r r ) is given. Going upstream, that is to 
the left in the diagram in analogy with equations (155) and (156), the desired velocity VF^{v r ) and the 
desired acceleration VF^(a r ) are obtained. The force generator Fff r is analogous to /2 in figure 33, 
except that Fff r may not be of the pure feedback type. To proceed further to the left, it is necessary 
to obtain a relative inverse of Fff r . But relative to what? That is defined by a new discrete variable, 
namely the control mode, (i c • For example, suppose that the test aircraft, unlike the Harrier, has in 
addition to throttle and nozzle a side force controller, so that the power control u p has complete, 3-D 
control of force. Then the aircraft rotation and translation may be decoupled within reasonable limit. 
That is, in addition to specifying the motion of the center of mass by VF 4 (r r ), the desired rotation 
can also be specified. Then, assuming that the configuration VF^(uf) and other parameters 
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VF 2 (p ) are given, the inverse of Fff r relative to VF 2 (u p ) and the inverse of the moment generator 
F™ 6 relative to the moment control VF 2 (u m ) can be generated. In particular, the controls ( u m ,u p ) 
are given by the six-dimensional static inverse of both moment and force functions. That is, 


Ufn Ffr {jTl b , df , Uj , Tf , ty , Cfyp , ^brbi P) 

Up = Fb r {™ b , a r, u },r r ,v r ,C iyr ,u brb ,p) 

so that 

rn b = F^ium, u p ,Uf,r r ,v r , 
a r — F b r (u r ni‘u p ,uj,rr,v r ,C br ,uJ brb ,p) 


(191) 


(192) 


Inverses such as those in equation (191) are sometimes referred to as trim maps. This particular zero- 
order trim map can be raised to any order by means of dynamic forms, as discussed previously in the 
algorithm (149). At this point the dynamic inversion of the model in figure 35 has been completed by 
obtaining VF 2 (u m ) and VF 2 (u p ). The inverse model is shown in figure 37. This control mode will 
be designated by the mnemonic p c = DFC for “direct force control.” 



Figure 37. Direct force control (DFC) mode. 


Another very useful control mode is the direction cosine control, p c = DCC. For example, the first 
step in the recovery from a severe atmospheric disturbance may be to recover the aircraft attitude and, 
only after the attitude stabilizes, to proceed to recover the trajectory. In this mode the desired rotation is 
specified at each instant by the rotation form RF 2 (C br ); the power and configuration controls VF 2 (u p ) 
and VF 2 {uf ) are driven according to a preassigned schedule. In this mode the trajectory VF 2 (r r ) is 
a dependent variable. The corresponding trim map is simply the 3-D moment trim: 

Um = F b lb (m b ,u p ,Uf,r r ,v r ,C br ,u) brb ,p ) 

so that 

m b = F^ b (um,u p) Uf,r r , Vr, Chr.oj^p) 
a r = F 6 or (u Tn ,'i/p,u^,r r ,v r ,C' fcr ,a; frr6 ,p) 

The composite algorithm producing force subject to moment trim will be denoted with an overbar, as 
follows: 

a r = F b r (m b , u p , u f , r r , v r ,C hr ,u brb ,p) = F 6 ar o F 6 m6 (195) 


(193) 

(194) 
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where the last term means that the algorithm F™ b is called first and then algorithm F^ r . The DCC 
control mode is useful, not only for attitude maneuvers but also for the design of other guidance laws. 
This mode may be interpreted as a nonsingular change of control coordinates from the natural moment 
control variable u m to the new control variable namely, angular acceleration. The corresponding 
aircraft model is shown in figure 38. The natural moment control u m is a dependent variable, which is 



Figure 38. Direction cosine attitude control (DCC) mode. 


computed by the moment trim algorithm F™° . The DCC mode provides a very convenient and realistic 
procedure for the suppression of zero dynamics. A pure feedback approximation of the force model 
may be obtained by evaluating the moment-trimmed force with m 5 = = 0 , namely, 

a? = F§ r (Q,u p ,Uf,r r ,v r ,C br ,Q,p) (196) 

The error due to the suppression of the effects on a r of u m and ui ^ is 

^or = F\) Up, Uj , Tf , Vr > Cfj r , ujfoj-h, p) Ffo (0^ Up, Uj , Vi', Vr , C^ r , 0, p) (197) 

The present report is concerned with cases for which this error is small. 

Still another very useful mode is based on the alpha control mode shown in figure 39. Consider 
the Euler angle control (EAC) mode p, c = FAC , in which the moment function is inverted with 
respect to the moment control u m . 

u m = Ft na (rn a ,Up,Uf,r r ,v r ,att, t ,ab t ,Ctr,utrt,utrt,p) (198) 

The result is shown in figure 39. The composite algorithm producing the moment-trimmed force is 
denoted, as before, with an overbar 

pm a = F ft o pma (199) 
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Figure 39. EAC mode. 


The pure feedback approximation of the force is 

ft = Ft na (°» u f> r r,v r ,oc bt , 0 , Ctr , 0 , 0 ,p) ( 200 ) 


The action of the control mode may be described in general as follows. Consider the expanded 
force and moment algorithm Ft, as defined in equation (188), with input xt and output yt. 


Vt = Ft(x t ) ( 201 ) 

Suppose that some coordinates of yt are labeled as independent ‘i’ and the remainder as dependent ‘d.’ 
Similarly, let some coordinates of xt be labeled as active controls ‘c’ and the remainder as parameters 
‘p.’ Each such labeling defines a control mode, provided that the number of independent coordinates 
equals the number of active control coordinates, which is the dimension of the control mode. The 
algorithm Ft for computing the corresponding partial inverses is given by 

(x t ,y t ) = F t {x * t , yt , y c ) 

M-F t (x t )L c < e ( ^ 

The first input x* contains the initial values for the coordinates of xt that are designated as active controls; 
the remainder of xt contains the parameters p. The second input yl contains the target (desired) values 
for the independent coordinates of yt, the remainder of yt containing the dependent coordinates is 
ignored. The third input pt c , the control mode, designates which coordinates are which. The algorithm 
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fills in the locations of xt designated as controls and locations of yt designated as dependent. Table 4 
contains several practical examples of control modes. 


Table 4. Example control modes 


Ac 

m b 

m a 

m r 

ft 

h 

Clj" 


u m 

Up 

u f 

a 

SIM 

ddd 

ddd 

ddd 

ddd 

ddd 

ddd 

ddd 

ppp 

ppp 

ppp 

ppp 

DCC 

Hi 

ddd 

ddd 

ddd 

ddd 

ddd 

ddd 

ccc 

ppp 

ppp 

ppp 

EAC 

ddd 

Hi 

ddd 

ddd 

ddd 

ddd 

ddd 

ccc 

ppp 

ppp 

ppp 

DFC 

Hi 

ddd 

ddd 

ddd 

ddd 

Hi 

ddd 

ccc 

ccc 

ppp 

ppp 

HAR 

ddd 

Hi 

ddd 

Hi 

ddd 

ddd 

ddd 

ccc 

cpc 

ppp 

ppc 

VAT 

ddd 

Hi 

ddd 

Hi 

ddd 

ddd 

ddd 

ccc 

cpp 

ppp 

ccp 

HEL 

ddd 

Hi 

ddd 

Hi 

ddd 

ddd 

ddd 

ccc 

cpp 

ppp 

cpc 

CFY 

ddd 

Hi 

ddd 

Hi 

did 

ddd 

ddd 

ccc 

cpp 

ppp 

ccc 

CFB 

ddd 

Hi 

ddd 

Hi 

ddd 

ddd 

did 

ccc 

cpp 

ppp 

ccc 

CFF 

ddd 

Hi 

ddd 

Hi 

did 

ddd 

ddd 

ccc 

cpp 

cpp 

pcc 

CFX 

ddd 

Hi 

ddd 

iid 

did 

ddd 

ddd 

ccc 

ppp 

ppp 

ccc 

CFZ 

ddd 

Hi 

ddd 

dii 

did 

ddd 

ddd 

ccc 

ppp 

ppp 

ccc 


Example 1: p c = SIM, the simulation mode. In this mode all coordinates of input xt are 
parameters, and all the output yt is dependent. There is nothing to invert. That is, SIM simply invokes 
the forward map Ft. 

Example 2: p c = DCC , the previously described direction cosine attitude control mode. In this 
mode, m b is designated as independent and u m is the active control. 

Example 3: p c = EAC, the previously described Euler angle attitude control mode. In this mode, 
the independent output is m a and the active control is u m . 

Example 4: p c = DFC, direct force control mode. As already described, this mode requires 
a six-dimensional inversion. Both angular and translational acceleration (mj,a r ) are given, and the 
inverter finds the required (trim) values of the moment and power controls ( u m , u p ). 

If these four modes were the only ones of interest, the trim problem would be solved. But, in fact, 
most control modes of practical interest lead to more-complicated inversion problems. Several such 
cases are described. 

Example 5: p, c = HAR. This is the Harrier mode, in which the direct side force in not available. 
It is treated as a parameter, u p 2 = 0. The active control is five-dimensional: full moment control, 
plus throttle and nozzle. The active control has direct access to first and third components of the body 
coordinates of the total force, f b . Full control of f r is achieved by tipping the aircraft. That is, one of 
the attitude degrees of freedom is used for control. In the HAR mode, that degree of freedom is the 
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roll (t*to 3 ) channel, while the pitch (a bt i) and yaw (a bt2 ) channels remain free parameters. The relative 
attitude is represented by C bt = E 2 (a bt i)E^(a bt2 )Ei(a bt3 ). 

Example 6: /x c = VAT. This is the vertical attitude control mode, in which the aircraft is sitting 
on its tail. The force is controlled by the thrust and the pitch and yaw channels. The force balance must 
be achieved by means of the throttle u p \ and the relative pitch and yaw angles a bt i and a bt2 , while the 
relative roll angle a bt 3 acts as a parameter: in the VAT control mode the aircraft may execute an inde- 
pendent roll schedule c* 6 t 3 - The relative attitude is represented by C bt = E 2 (a bt i)E^(a bt2 )Ei(a bt ^). 

Example 7: n c = HEL. This is the helicopter mode, in which neither the direct side force 
controller nor the nozzle is available for active control. Only the full moment control and throttle 
(collective in helicopters) are active. Two additional degrees of freedom for complete force control are 
provided by roll and pitch channels, while the yaw channel remains free. The force balance is achieved 
by manipulating the throttle (or collective) u p \ and the relative pitch and roll angles a bt i and a bt 3 , 
respectively, while the relative yaw angle a bt2 is treated as a given parameter. 

Example 8: n c = CFY . This is the conventional (high-speed) cruise control mode. The thrust 
is mainly in the direction of the longitudinal axis b\ of the aircraft. The direct force control is the 
throttle. All three attitude channels are used for control purposes: pitch and roll to provide the two 
additional degrees of freedom for complete force control. Yaw is used to maintain the body side force 
f b 2 near zero. The seven-dimension inversion computes the moment control, throttle control, and the 
three relative Euler angles: (u m ,u p i,a bt ). 

Example 9: The CFB mode is also seven-dimensional, but it is the side air velocity v b2 that is 
externally specified (usually zero). 

Example 10: In the CFF mode, flap ut\ is manipulated to obtain a given angle of attack a bt i- 

The last two entries in the table are modes for which the throttle is not available for control; it is 
being moved between its limits or sitting on one of them. 

Example 11: The CFX, conventional flight X force control mode, is essentially an airspeed 
control mode. The independent force components are 

( ftl > St2i fb2 ) ( 203 ) 

The normal force channel fa is not controlled; it is a dependent variable so that the climb rate h is not 
controlled directly. 

Example 12: On the other hand, in the CFZ mode the independent force components are 

( ft2i fb2 ) (204) 

while the longitudinal channel ft\ is not controlled. This is a climb rate control mode in which a given 
altitude rate h(t) schedule is being executed, while the airspeed is not directly controlled. 

These examples illustrate the fact that there may be many control modes to consider in a complete 
design of the guidance subsystem. Each mode generates a specific inversion problem. Furthermore, it 
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is desirable to be able to easily add new control modes, without having to redesign the whole algorithm. 
In the present approach, a new mode may be added by simply entering a new line in the control mode 
table. We use a standard variable dimension Newton-Raphson inversion algorithm. The labels in the 
control mode table are used to rearrange the input to this standard algorithm by storing the coordinates 
of labeled as independent in the head yj of the standard output, y*\ the dependent part is placed in 
the tail y* D of y* . Similarly, the control coordinates of x * are placed in the head x* c of the standard 
input x*, and the parameters in the tail x* p . For a valid control mode, the Jacobian dyj/dxc is always 
nonsingular in the region of interest. When the standard inverter reaches the threshold, then the output 
(x,y) is redistributed according to the labels into ( xt,yt ). If the partition functions are defined by 

(x c ,x P ) = n x (x t ,p c ), x t = iv~ l (xc,xp,yL c ) 

{vhVD) = ^i i{xt,vc), yt = ^y l {yi,yD^c) 

Then the multimode (direct) function for the standard inverter is given by 

F{x,Hc) = 7 Ty { Ft [7T~V,/Z C )] ,// c } 

In this fashion a single static inversion algorithm for any valid control mode is obtained, 

{xt, yt) = F t {x* t , yt,Hc) (207) 

The parameter e is introduced in order to be able to suppress zero dynamics, if any, as follows. 


(205) 


(206) 


4(4.S/?.#*,<0 


Ft{x$,yt,Hc) if e = 1 

F t {x* t , yl , Mc)m b =u brb =0 if f = 0 


(208) 


That is, F} (x$ , y* , jj, c , 0) computes the force subject to the constraint that both angular acceleration t 
and angular velocity c of the aircraft relative to the runway are zero. 


The Ft algorithm is quite powerful. For example, 

{xt, yt) = f}{x*, yt, SIM, 1) (209) 

reproduces the exact simulation. In this trivial case, there is only one call to Ft and no inversion. The 
input x* is copied without change into xt. The input y% is ignored, while the output yt = Ft{xt). 


The pure feedback approximation of the forward map Ft corresponding to equation (200) may be 
obtained with two calls: 

{x t ,yt) = F}{x* t ,y* t ,SIM, 1) 

(x?,y?) = F}{x t ,y t ,EAC, 0) 

The m a component of yt is the exact second derivative of the Euler angles, m a , while the ft component 
of y^ is the pure feedback approximation / ( °. The resulting force error 

ft = ft -ft (211) 

In the absence of zero dynamics, that error is bounded by the trim error bound e. This inverter error 
will henceforth be largely ignored. 
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(212) 


The pure feedback inverse of the force generator is obtained with two calls: 

(*?■»?) = Fl(x;, y ;,vc,o) 

t 0 t * 

x\ = m a = m a 

(it, Si) = F}(x\, V ' t ,EACfi) 

The first call trims the aircraft to the desired force / t ° = / t * by manipulating the active force controls 
as specified by the control mode n, but with angular acceleration and angular velocity held at zero. 
The resulting moment control is u® n , and the power and configuration controls are (u°,u^). Then x ^ 

is copied into x\, except that x\ contains the true angular velocity d^. The second call trims the 
moment control for the specified angular acceleration ra* and angular velocity d£ f , disregarding force 
balance. The overall result is the exact moment control u m for the given angular acceleration and 
angular velocity, the pure feedback values u°j, and for the given / t *, and the exact f t . The force 
error due to the neglected zero dynamics is, again, 

ft = ft ~ft (213) 

assuming a perfect inverter. 

The final step in the construction of the pure feedback approximation of the multimode dynamic 
inverse T° and direct T° is to raise the (zero-order) f} algorithm to order two. As described previously, 
that step may be done routinely by means of dynamic forms. The result is shown in the following block 
digrams. Figure 40 shows the pure feedback approximation T° of the exact inverse T. The function 



pure feedback exact 


Figure 40. Pure feedback approximation T° of the inverse map T. 
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F t is the multimode static inverter, as defined by equation (208) raised to arbitrary order. The input to 
the inverse map is the mode //, e = 0, which cuts the zero dynamics, order k — 2, and second-order 
dynamic form of x*, namely DF 2 (x*), and second-order vector form VF 2 (y*). The output of this 
call to Ft is the filled-in DF 2 (x°) and VF 2 (y°). The next step is to copy the zero-order components 
DF°(x°) into £>F 0 (xj), and copy the first and second derivatives of a into the appropriate slots of 
x\, and similarly with u t rt and u tr t- Finally, the last call to Ft performs the exact moment trim. At 
this point the pure feedback inverse T° has been constructed, since the given dynamic forms have been 
transformed into the corresponding state and control. 

The pure feedback direct algorithm T°, shown in figure 41, is the reversal of the preceding algo- 
rithm. The first step is to compute the exact angular acceleration from the given state and control. 
Next, the second-order angle form AF 2 (c* 6t ) is assembled, and the pure feedback static inverter is 
called first with order zero. The resulting acceleration a r is loaded into VF 2 (r r ) and VF l (v r ). Then 
Ft is called again, but to order k = 1, and the process is repeated once more. Finally, the output, 
DF 2 (x £*) and VF 2 (y t °), is computed. At this point T° has been constructed, since the state and control 
have been transformed into the corresponding dynamic forms. Note that the approximate transformation 
T° suppresses (blocks) the effects of zero dynamics. The effects of these two transformations may be 
summarized with the help of the block diagram in figure 42. 



exact pure feedback 


Figure 41. Pure feedback approximation T° of the direct map T. 
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Figure 42. Pure feedback approximation for aircraft guidance. 


The system 

x = f(x,u,p, 0) (214) 

is viewed through T° and T° as a decoupled string of integrators (see bottom of figure 33). There will 
be deviations from this Brunowsky form if there are zero dynamics present. The control of such errors 
is considered next. 

4.5 Guidance Servo 

As discussed previously (see figure 32), a simple and direct way to control the errors caused by the 
neglected zero dynamics is by means of the servo structure shown in figure 43. The right-hand block 





guidance regulator 


Figure 43. Guidance servo in canonical coordinates. 
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represents the exact simulation, whose state and control (x 9 , u 9 ) are taken as the solution to our guidance 
problem. By construction, (x 9 ,u g ) always satisfies the state equation. The tracking error is given by 
the difference between the commanded output produced by the guidance output trajectory generator, as 
shown in figure 5 and discussed in Chapters 1 and 3 , and the actual output of the simulation: 

VF 3 (yf)VF l (y%) (215) 

For the pure feedback case, these errors are identically zero. In any case, the loop on that error is 
closed by assuming that the system being regulated is a pure Brunowsky form. The corrective feedback 
is injected into the fourth derivative of y\ and the second derivative of y%, which are the canonical 
controls for the pure feedback case: 


(yf)( 4 ) = K?(yf,yf,yf,(yf)^) 

y e 2 = Km,y$) 


(216) 


The resulting modified signals VF A (y x ) and VF 2 (y 2 ) are then transformed as before by the pure 
feedback dynamic inverse f° into the natural state and control. The state is used only to initialize the 
simulation while the control drives the simulation model. 


The design of the regulator law K 9 is greatly simplified by the fact that the system to be controlled, 
when viewed through the transformations T° and T°, is, except for the corrupting influence of the zero 
dynamics, a collection of simple, decoupled strings of integrators. Linear techniques such as pole 
placement are quite effective. Of course, the regulator may be nonlinear or include dynamic elements 
such as filters to obtain smoothing. No gain scheduling is needed. That function is being provided 
automatically by the linearizing transformations. 

In the following numerical examples, a linear regulator was selected and all the gains were held 
constant throughout the maneuvers. For each fourth-order string of integrators, the closed-loop poles, 
including a first-order filter with time constant of r = 0.18, were placed at 

(wi.Cl) = (0.5, 0.7) 

(W2.C2) = (2.0, 0.7) (217) 

a = 2.0 

For each second-order string of integrators, the closed-loop poles were placed at 

(wi,Cl) = (2.0, 0.7) (218) 


and filters were not used. 
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5 EXAMPLES 


In this section we present two numerical examples. 

5.1 Maneuver: Move Forward 300 Feet 

Consider again the maneuver M\ described in figure 10. The commanded motion of the center of 
mass, VF 4 (r£), which transfers the Harrier 300 feet forward in 18 seconds, from hover to hover, is 
shown in figure 44. There are four segments. 1 ) In the first segment of the maneuver, 0 < t < 8, the 
aircraft is accelerated to 30 ft/sec in the Harrier mode, namely n = HAR, while the pitch angle, which 
is a parameter in this mode, is held at zero degrees. 2) In the second segment, 8 < t < 10, the aircraft 
coasts for two seconds at 30 ft/sec with n = HAR. 3) In the third segment, 10 < t < 18, the aircraft 
decelerates to hover with /j, = HEL, while the nozzle is commanded to 90 degrees. In this mode the 
pitch axis is used for control, while the nozzle angle becomes a parameter. 4) In the last segment, 
18 < t < 24, the command is to hover in the HEL mode. 



Figure 44. Forward step: commanded position VF 4 (r c ). 


The behavior of the open-loop controls is shown in figure 45. The first panel shows the behavior 
of the pitch moment control u m 2 . The solid curve is 2 , which is the position of the pitch moment 
control assuming zero angular velocity and zero angular acceleration. The dotted curve is u m 2 , which 
is the position of the pitch moment control for the given angular velocity and acceleration without 
correcting for the effects of zero dynamics on the resulting force. The second panel shows the throttle 
and two of its derivatives, VF 2 (u p 1 ). The third panel shows the nozzle angle and two of its derivatives, 
VF 2 (u p 3 ). The last panel shows the motion of the relative pitch angle and its derivatives, VF 2 (otf, t i). 
This is also the total pitch angle 6 of the aircraft body axes since in this case the t-axes were chosen to 
coincide with the runway axes, namely, Ct r = I- 
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0 5 10 15 20 25 

time (sec) 


Figure 45. Forward step: open-loop controls. 


The t-coordinates of the force error caused by the zero dynamics is the difference 

ft = ft~ ft (219) 

fr 

between the exact force ft, computed with the exact (e = 1) force algorithm F( , and the commanded 
force 

ft=Ctr(yi/g-h) (220) 

This force error, shown in the first panel of figure 46, is quite small. The small jagged sections are 
inverter errors. The rest is due to zero dynamics. Some cross coupling from the longitudinal to normal 
channels is evident in the HEL mode. The RCS jets that generate the pitch moment also generate normal 
force. The inverse map by design does not take this effect into account. The resulting translational 
acceleration errors are bounded by 0.25 ft/sec. 

Thus far in the present example the open-loop behavior has been considered. Now consider the 
closed-loop behavior, as in figure 43. The closed-loop tracking error rf is shown in the second panel 
of figure 46. The maximum error for the entire maneuver is only 0.2 feet along track, and 0.6 feet 
in altitude. The acceleration error a® = rf, shown in the third panel, is bounded by the small value 
of 0.25 ft/sec 2 . The corrective feedback, (yf)( 4 ) = (rf )( 4 ( (see figure 43), shown in the fourth panel, 
is bounded by the small value 0.25 ft/sec 4 . So, the weak zero effects, in the present case of the step 
forward maneuver, are well controlled by the simple and rather crude (guidance) regulator. 
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Figure 46. Forward step: force errors. 


5.2 Maneuver: Frisbee 

As a second example, consider the Frisbee maneuver described in figures 12 and 13. This maneuver 
is synthesized by means of 14 segments. During this complex maneuver, both the control and the attitude 
coordinatization modes change. As shown in the first panel of figure 47, the trajectory transfers the 
center of mass from hover at the origin to 4000 feet along the runway and 667 feet in altitude. The 
second panel shows the behavior of the velocity v£. The duration of the maneuver is 50 seconds. The 
horizontal velocity changes from zero to 215 ft/sec, and the vertical from 0 to 15 ft/sec climb rate. 
The lateral velocity remains at zero. The third panel shows the yaw angle. Note the 360-degree yaw 

maneuver. The large-angle maneuver is accomplished by means of an intermediate axis system t , as 
described before. In the first part of the maneuver, Ctr is parameterized by means of Euler angles in 
the standard sequence: 

C tr = Ei(<fi)E 2 (6)E 3 (il>) (221) 

Then during the second, higher speed part, the wind path axis (WPA) parameterization is used. 

The last panel shows the behavior of the body coordinates of the relative air velocity v% c . It may 
be noted that during the initial, slow part of the maneuver, there is significant side velocity, and, in fact, 
for a time the aircraft is flying backward, but in the high-speed condition the speed is along t\ with no 
side velocity. 

The behavior of the open-loop controls is shown in figure 48. The first panel shows the moment 
controls u m . There are three pairs of curves. The difference between members is due to the zero 
dynamics. The next two panels show the throttle and nozzle and two of their derivatives. The last panel 
shows the flap and aileron droop. 
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The open-loop behavior of the angular variables is shown in figure 49. The first part (the first 
20 seconds) of the maneuver is accomplished in the HAR mode, and the rest in the HEL mode. 

4 — ► 

The first panel shows the motion of the relative angles a i t locating the body axes b relative to the 

intermediate axes t . In the first part of the maneuver, the roll angle a ^ rolls the force vector to 
maintain along-track acceleration. In the second part, the relative pitch angle a^i tips the force vector 
forward to maintain along-track acceleration. Later in the maneuver, this angle becomes the angle of 
attack (producing lift) because the coordinatization of aircraft attitude is switched to the wind path axes. 



Figure 49. Frisbee: open-loop rotation and force error. 


The second and third panels in figure 49 show the first and second derivatives of a^ t . Both vectors 
are well behaved. The force error ff [see eq. (219)], due to zero dynamics, is shown in the bottom 
panel. It is quite small, being bounded by 0.005 g’s. The activity in the HAR mode is due to the RCS, 
whereas the activity in the HEL mode is due to both RCS and the lift of the elevator. 

The closed-loop position error rf, shown in the first panel of figure 50, is bounded by 1 ft throughout 
this complex maneuver. The acceleration error a®, shown in the second panel, is quite well behaved, 
being bounded by the small value of 0.3 ft/sec 2 . The corrective feedback (yf)( 4 ) (see fig. 43) shown 
in the third panel is active, but bounded by the small value of 1 ft/sec 4 . The final panel shows the 
closed-loop behavior of the moment control u m . 
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Figure 50. Frisbee: closed-loop behavior. 


In summary, the two fairly complex examples show that the proposed approach based on inversion 
of the pure feedback approximation, followed by a regulator in canonical coordinates, is feasible for 
maneuvers that are slow relative to the aircraft zero dynamics. 
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6 CONCLUSION 


An effective method has been described for guiding a dynamic system, such as an aircraft, through 
complicated maneuvers. The maneuvers are specified by means of a control point table, which is the 
discrete input to the guidance algorithm. The output is the continuous reference motion of the aircraft 
state and all the controls. The state is maintained continuous across the control points. Furthermore, 
this reference motion is flyable exactly by the aircraft simulation model. Therefore, the structure of 
the regulator for the actual aircraft may be taken to be that of an exact model follower. This feature 
has the great practical advantage that this (plant) regulator is subjected to process uncertainties only 
and is not loaded down by apriori known signals. Consequently, the full authority of the regulator 
is available for the control of disturbances and other uncertainties. The proposed multimode dynamic 
inversion procedure, consisting of the pure feedback approximation followed by a (guidance) regulator, is 
effective for maneuvers that are slow relative to the the plant zero dynamics. For stronger zero dynamics, 
the method described in reference 10 for the improvement of the pure feedback approximation may be 
employed. The methodology of the dynamic forms is very effective for organizing and drastically 
simplifying the coding of the proposed algorithms. 
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